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Methods of Applied Dynamics 


Description: 

This monograph is designed to give the practicing engineer a clear understanding 
of the principles of dynamics with special emphasis on their applications. Begin- 
ning with the basic concepts of kinematics and dynamics the course proceeds to the 
discussion of the dynamics of a system of particles. The analytical (Lagrangian) 
method of dynamic analysis is treated in full detail. Both classical and modern 
formulations of the Lagrange equations including constraints are discussed and ap- 
plied to the dynamic modeling of aerospace structures using the modal synthesis 
technique. A list of references is given at the end of the monograph. 
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Chapter 1 
Kinematics 


Kinematics relates to the geometry of motion disregarding the forces causing the 
motion. 


1.1 Vectors 

A vector has direction and magnitude (velocity, force, etc.) 

Physical types of vectors : 

1. free vector: velocity 

2. sliding vector: force on rigid body 

3. bound vector: position, force on elastic body 

NOTE : All mathematical operations with vectors involve only their free vector 
properties of magnitude and direction. 

Addition : The sum of two vectors is represented by the diagonal of a parallelo- 
gram formed by the two vector sides. 
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NOTE : Vectors are denoted by bolding the symbol. 

Vector addition is commutative and associative. 

Unit Vector : A unit vector has unit magnitude (length). 

Vectors are often conveniently expressed in terms of unit vectors: 

A = -4iei + ^262 + A 3 e3 

where A?, A 3 are known as (scalar) components of the vector A. 
Often unit vectors are used which are (mutually) orthogonal : i,j, k 



A — . l x i + A y j + A z k 
Here the components A x ,A y ,A z 
are the orthogonal projections of A 
onto the coordinate axes x, y, z. 


NOTE: The unit vector in the direction of the vector A is identified as: 

e A = A/A A = |A| 

P-S. #1 : The German word for unit is “EINHEIT.” This is the origin of the common 
use of the letter e for the unit vector. 

P S. #2 : It is important to use suggestive symbols for physical quantities: m = small 
mass; M = large mass. 

Scalar (“Dot”) Product 

A • B = A B cos 0 

6 


1 



where 6 is the (smaller) angle between A and B. Expressed in terms of orthogonal 
unit vectors e 1; e 2 , 63: 


A • B = (Aj ej + A 2 e 2 + A3 63) • ( Bi ei + B 2 e 2 + B 3 03) 


— A; By A A 2 B 2 + A 3 B 3 


NOTE: 


e x • ei = e 2 • e 2 = e 3 ■ e 3 = 1 (UNIT LENGTH) 

d ■ e 2 = e t ■ e 3 = e 2 • e 3 = 0 (ORTHOGONALITY) 

Vector (“Cross”) Product 


A x B = A B sin 6 N 


O<0<7T 


where 8 is the smaller angle between A and B. N is perpendicular to A and B 
such that A, B. N form a right-handed system. (Right hand thumb rule) 

Expressed in terms of orthogonal unit vectors e 1 ,e 2 ,e 3 . 

A x B = (Ai ei + A 2 e 2 + A 3 e 3 ) x (B v e L +fi 2 e 2 + B 3 e 3 ) 

= 0i (A 2 B 3 — A3 fi 2 ) + 02 (A3 B\ — A 1 B 3 ) 

+ e 3(A] R 2 — A 2 B\) 

NOTE: 


ei 

x e! 

= 0 2 X 0 2 

= 0 3 X 0 3 

ei 

x e 2 

= C 2 x 0 ! 

= e 3 

e 2 

x e 3 

= -0 3 X 0 2 

= ei 

e 3 

x ei 

= -e, x 03 

= e 2 
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Mnemonic Code: 


e 3 



Arrange the three unit vectors ei,e 2 ,e 3 
in clockwise order. The vector product of two unit 
vectors is equal to the third if the two follow each 
other clockwise and equal to the negative third if 
the two follow each other counterclockwise. This rule 
is especially useful if some components are zero; e.g. 


a x b — (o[ e! + a 2 c 2 T 03 e 3 ) x fti ej 
= — a 2 bi e 3 + a 3 bi e 2 


NOT E: 


A x B = — B x A 


(anticommutative) 


Scalar Triple Product: 


A • (B x C) = B ■ (C x A) = C • (A x B) 


This product is geometrically equal to the volume of a parallelepiped of sides 

A, B, C. 

Vector Triple Product: 
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A x (B x C) = (A • C)B - (A • B)C 


Moment of a Vector: 

If the first vector in the vector product A x B represents a position vector R then 
the resultant vector product is called the moment of the vector B: 


M = (R x B) 


Velocity: It is the time derivative of the position vector R: 

V = A = R |V| = Speed 

|R| = Distance 

Acceleration: It is the time derivative of the velocity V: 



<f 2 R 

dt* 


= V = R 


NOTE: 

Sometimes we need the time derivatives of the scalar and vector products: 


H A • B) 

II 

> 

+ 

dA r> 

It U 


= A B 

+ 

A B 

U A x B) 

— a x — 

A x dt 

+ 

*b »B 


= A x B 

+ 

A x B 
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R = R(0 V = V(0 


1.2 Angular Velocity 

Finite angular rotations are not commutative and therefore cannot be treated as vec- 
tors. (Rotate a book 90° about the x and y axes and repeat the procedure in 
reverse order and observe the difference in final orientation). Angular velocity can 
be shown to be a vector. Consider the rotation of a point P about an axis 00 ; 
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We now assign a vector u> in the direction of 00' with the magnitude equal to uj. 
Then the linear velocity v of the point P due to this rotation is 

v = wxR R = CONST (1.1) 

The magnitude of v is clearly uj R sin 8 and its direction is normal to the plane 
spanned by uj and R. (Right-hand thumb rule: fingers indicate rotational sense, 
thumb is in direction of uj). If now a second axis is given through O , then the linear 
velocity to this rotational rate 0 is given by V = fi x R. The total velocity of point 
P is then 


V p = v + V = (wxR) + (/? x R) = (u + S?) x R 
= f2 p x R where fl p = uj + fi Q. E. D. 


1.3 Vector Derivative in a Rotating Frame 

A vector A is seen by an observer in a fixed frame A , Y\ Z and also by another observer 
in a rotating frame x,y,z. Unit vectors in the fixed frame are denoted by i. j. k 
and in the rotating frame by e l5 e 2 ,e 3 . The rotating frame has angular velocity oj 
relative to the fixed frame. 


A — A\&i + A262 + ^363 
= A x i + A y j + A z k 
A = generic vector 
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The rate of change of A as seen from the fixed frame: 


A — -4ie! + T 2 e 2 + A 3 e 3 + + A 2 e 2 + A 3 e 3 


But: 


ei — u> x ei , e 2 = w x e 2 , e 3 = w x e 3 


Therefore: 


A = (A), e / + w x A where (A), e( = Aie[ + A 2 e 2 + A 3 e 3 


( 1 . 2 ) 


This relation holds for any two systems ^4 and B : 


( A) a — (A )b + &ba x A 


NO I L - If A — ^ ba then — {wba)b where u*ba is the angular velocity of 

B as seen in .4. 


1.4 General Motion in a Moving Frame 


The X, Y, Z frame is fixed (inertial frame) and the x,y, z frame rotates relative to it. 
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Ro = Position of Origin 




R = R 0 +wxr + wx(wxr) + 2(wxv)+a 


(1.4) 


Alternate Form: 


R = V 0 + wxV 0 + «xr + wx(wxr) + 2(wxv) + a 


(1.5) 


NOTE : 

It is very important to have a thorough understanding of the physical meaning of 
the five acceleration terms on the right-hand side of Eq. (1.4). 

1. R 0 = acceleration of 0 of moving frame (D’Alembert/Einstein acceleration) 

2. w x r = “slingshot” acceleration (Euler acceleration) 

3. w x (w x r) = centripetal acceleration 

4. 2 (u7 x v) = Coriolis acceleration 

5. a = relative acceleration as seen in moving frame 

P. S. 

The Coriolis (1792-1843) acceleration is somewhat difficult to visualize. It is com- 
posed of two separate kinematical effects: one velocity change is due to a change 
in the direction of v due to (“slingshot” effect), the other velocity change is due 
to a radial change of the point P position. Both changes of velocity are equal to 
(w x v) resulting in the factor 2 in the Coriolis acceleration. Mathematically: 

a! = (u> x v) 2nd term of Equation 1.2 

a 2 = W X ( ~~ ) Rel = (W X V) 
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Tangential and Normal Components 


Consider the position of a point P as it moves along a curved path in space. 


Z 



s = distance along curve 


Velocity : 


where 



dr 

ds 


ds 


Tt =v * T 


dr 

ej = — = tangent unit vector 


Acceleration : 


dej ds . ,rfer 

a = v = ve T + v e T = ve T + v—— ■ — = ve T + v 

ds dt ds 


Define: 


dej 

ds 


= /ce A - 


( 1 . 6 ) 


15 



where 

k = curvature [rad/meter] 

- = p — radius of curvature 

n r 

e N = normal unit vector 
The acceleration is then given by 


a — a 6j -f* v /c6yv — o Gj — e \ 

P 


The first term is the tangential acceleration and the second term is the normal or 
centripetal acceleration. 

Define a third unit vector to complete the orthogonal triad of unit vectors (Trihe- 
dron) at the point P: 


e B = e T x e iV 


e B = binormal unit vector 


4 

5 

U 

0 

fl 

£ 


G B Normal plane 
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(er,ejv) = Osculating Plane 
(eB,eAr) = Normal Plane 
(er,e B ) = Rectifying Plane 

The curvature k measures the rotation rate of the normal plane as the point P 
moves along the curve. 

Now we differentiate the binormal unit vector with respect to the distance along 
the curve: 


de B 

ds 


= e' B = e' T x e N + e T * e^, = (ne N ) x e N + e T x e' N 


e 7 x e 7v 


Therefore, X e T . 


Since e B is a unit vector, we also have X e B . Therefore, must be parallel to 
Cjv- 


Define : 


deg 

ds 


-t e N 


where r = torsion 


[rad/meter], \ = <t = radius of torsion 


(1.7) 


A positive torsion (r > 0) corresponds to a clockwise rotation for a motion of P 
along the curve. The torsion r = 0 for a plane curve. The torsion measures the 
rotation rate of the osculating plane as the point P moves along the curve. 


Since e^r = x ej we obtain the spatial derivative of e^v as: 


e 7 v = e*B x e T + e B x e' T = (-re N ) x e T + e B x (/ce^) 


or 


t/6yv 

ds 


= res — /cej 


DARBOUX Vector 


( 1 . 8 ) 


The set of Equations 1.6, 1.7, and 1.8 is collectively called the Frenet-Serret 

formulas. 

The angular velocity of the trihedron can be easily obtained in terms of the curvature 
yc, the torsion r and the velocity v of the point P as: 
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w = /c ueg - t v ej 


The first term represents the angular velocity of the tangential unit vector and the 
second term the angular velocity of the binormal unit vector. 

Ballistics Equations 

The motion of a projectile through the air is often analyzed using path variables, 
which are measurements made along the tangential and normal direction of the 
trajectory or path. Because of their convenience the (N-T) coordinates are referred 
to as natural coordinates . 

The resistance (drag) is taken proportional to a power of the tangential velocity v 


R = k f(v) 



The equations of motion are then 


m v = — mg sin 9 + k f(v) 


(1.9) 


mv 2 

= g m cos 6 

P 


( 1 . 10 ) 
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The radius of curvature can be related to the path angle 9: 


1 _ M _ 1 d9_ 

p ds v dt 


(1.11) 


The minus sign is taken to agree with our definition of the flight path angle 9 in 
the figure. 

Dividing by m and introducing Equation 1.11 in Equation 1.10 yields: 


i ; = —g sin 9 — c f ( v ) 


( 1 . 12 ) 


vO — —g cos 9 


(1.13) 


where c = — = ballistic coefficient 

To obtain the position of the projectile as a function of time in the (x, y ) coordinate 
system we have to use the inertial velocities: 


x — v cos 9 


(1.14) 


y = v sin 6 


(1.15) 


Thus, Eqs. (1.12) - (1.15) represent the equations of motion of the mass center of 
the projectile in the plane. Their general solutions have to be obtained by numerical 
integration. 
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Chapter 2 

Dynamics of a Particle 

2.1 Newton’s Laws 

We introduce Newton’s laws of motions as axioms. (Axiom = a self-evident or accepted 
principle.) The truth of these axioms is established by experimental verification or 
prediction. 

Isaac Newton (1642-1727) published his laws in 1687 in Latin. Using modern 
language they are 

1. Every body continues in its state of rest or of uniform motion in a straight 
line, unless compelled to change that state by forces acting on it. 

2. The time rate of change of linear momentum of a body is proportional to the 
force acting upon it and occurs in the direction in which the force acts. 

3. To every action there is an equal and opposite reaction; that is, the mutual 
forces of two bodies acting upon each other are equal in magnitude and op- 
posite in direction. 

In a rigorous sense these laws apply only to a mass point or single particle. 

NOTE: 

The first law is only a special case of the second law when there are no external 
forces. The third law will later allow the transition from the dynamics of a single 
particle to the dynamics of a system of particles. 
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Dimensions and Units 


The world of dynamics can be described in terms of four dimensions: Length (L). 
Time (T), Force (F), and Mass (M). It is, however, customary to treat the dimension 
of mass as a primary dimension and derive the dimension of force via Newton’s law 
F = ma or vice versa. 

NOTE : This contrivance is also common in other areas: distance expressed in 
carhours or lightyears. (Find other examples). 

SI - System (Metric) 


Basic Dimension 

Unit 

Length 

Meter (m) 

Time 

Second (s) 

Mass 

Kilogram (kg) 

Derived Dimension 

Unit 

Force 

Newton (kg m 


Customary System (British) 

Basic Dimension 

Unit 

Length 

Foot (ft) 

Time 

Second (s) 

Force 

Pound (lb) 

Derived Dimension 

Unit 

Mass 

Slug (lbfts~ 2 ) 


NOTE : 
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Weight = Force 
W = mg 0 

g 0 = 9.81 ^ = 32.2 f -\ 

s 2 s 2 

Question : 

a) How much does 1 slug weigh? 

w = mg 0 = 1 lb . SeC " X 32.2 = 32.2 /6s 

// sec" 1 

1 slug = 32.2 lb 

b) How much does 1 kg weigh? 

m kq m 

W* = mg 0 = \kg x 9.81 — = 9.81 

s 1 s z 

1 kg == 9.81 Newton 

P. S. The metric system sometimes allows for the auxiliary unit of force called 
kilopond (kp). 

1 kp = 9.81 Newton 

The customary system is sometimes used with inches as the length unit. The unit 
of mass is then variously called SLINCH, SNAIL or MUG (lb sec 2 /in). 

1 SLINCH=386.41bs. 
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2.2 D’Alembert’s Principle (1747) 


With a stroke of genius d’Alembert (1717-1783) wrote Newton’s law in the form: 


F + (— m R) = 0 


It is exactly this apparent triviality which makes D’Alembert’s principle such an 
ingenious step. D’Alembert principle introduces a new force, the force of inertia and 
makes possible the use of moving reference frame. The inertial forces act exactly like 
all the other forces. They cannot be distinguished in their nature from any other 
impressed (external) force. If an observer is not aware that he is in an acceler- 
ated system, then purely mechanical observations cannot reveal to him that fact. 
Einstein revised D’Alembert’s principle to a general principle of nature in his grav- 
itational theory (Equivalence principle). 

Einstein s “box experiment” (Gedanken experiment) 

A closed elevator is pulled upwards with constant acceleration g 0 while at the same 
time, the gravity field disappears. It is then impossible to distinguish between the 
following two hypotheses: 

1. The elevator moves upward with constant acceleration g 0 . No gravity field 
exists. 

2. The elevator is at rest but there is a gravity field of magnitude g 0 . 

NDIE: 

The name “apparent force” for the force of inertia is misleading, if it is interpreted 
as a force which is not as “real” as any other external force. Sometimes inertia 
forces are also called “fictitious” or “effective” forces. 


Because of D’Alembert’s principle it is possible to analyze all dynamic phenomena 
in making reference frames with strict rigor and consistency. 

Staying consistently in a moving reference frame, Newton’s law can be reformulated 
using Equation 1.4; notice that the interest lies now in the relative acceleration a 
and not in the absolute acceleration R. 
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ma=F - m R 0 — m (u> x r) — m[u> x (w x r)] - 2 m(u> x v) 


( 2 . 1 ) 


or 


ma = F + A + E + C + K (2.2) 

where 


A = — m R 0 
E = — m (u> x r) 

C = -m [w x (u x r)] 
K = -2 m( u) x v) 


d’Alembert force 
Euler force 
Centrifugal force 

Coriolis force or compound centrifugal force 


P. S. #1 

Sometimes it is stated that D’Alembert’s principle transforms a problem in dy- 
namics to one in statics. The Coriolis force does no work and is therefore called a 
workless (“wattless”) force. 

P. S. #2 

D’Alembert’s Principle may be called the Equal Rights Amendment (ERA) of dy- 
namics because it declares the inertial forces to be equal to any other force. 


2.3 Work; Kinetic and Potential Energy 


The work done by a force F as it moves along a path from A to B is defined as the line 
integral: 


W = 



F • dr 


Work 


With Newton’s law: 


F = m r 
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we obtain: 


W = 



F 


f B 

dr— m r ■ 

Ja 


d r 


Remember : 


i 1 d . , , . - l.n 

iT=--(r-r)dt = -dv 


Therefore: 



m 


m r ■ dr = — 
2 





Introduce: 



Kinetic Energy 


vr = t b ~t a 


( 2 , 3 ) 


The increase in the kinetic energy of a particle is equal to the work done by the 
external force. 

P.S. 

It was assumed that the force field F = F(r) is only a function of the position r 
and not of time also. 

Potential Energy 

Force Field F = F(r) 

I. The line integral /jf F • dr is independent of path. 
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The following alternate statements are all equivalent to I: 
II. The contour integral vanishes. 




dx =>■ 



d r = 0 


III. F = —grad V where V = Potential Energy 

IV. Curl F = V X = 0 V = Nabla operator (Del-Operator) 

NOTE : 

In Electrodynamics: E = — grad V where V = electrical potential. 

In Fluid Dynamics: v = grad <f) where (}> = velocity potential. 

The work done by the external force is now: 


W = 



F 



gradV • d r = — 




dV = V A -V B 


Inserting in Equation 2.3 yields 


V A - V B = Tb-T a 


or 


I A + T A — V B + Tg 

This is the principle of mechanical energy conservation . Therefore such forces or 
force fields are called conservative. 
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NOTE : 


Sometimes forces in nature are derivable from a time-dependent potential V = 
V(r, t). For these statements I to IV hold if the time is kept constant (t = CONST). 
These forces or force fields are called irrotational . Energy is noi conserved. 

P.S. 

In earlier practice, it was customary to use the negative of the potential energy 
which was called the work function U = —l'. In view of the above described 
conservation law it was an advantage to change this sign. The operational “vector” 
V was introduced by Sir William Hamilton (1805-1865). The name “nabla” was 
coined by Oliver Heavyside (1850 - 1925) after an ancient Assyrian harp whose form 
it resembles. It is also called “atled.” (Delta spelled backwards) 

Consider now the work done against the external force: 

t B [B ( B 

W = - F • dr = / grad V ■ dr = dv = V B - } ' A = AV 


Therefore, the difference of the potential energy is the work I have to do against 
the force going from A to B. 

Example 1 : 

Gravitational Force 


F = - 


mg Rg 

R 2 


e r R 0 = earth radius 


Work done against gravity: 

W = ~ j F ’ dT = m 9 *o / = rngRH- - -) = V(r) - V(r 0 ) 

Jr o it Tr\ T 


At the reference point r 0 we set the potential energy to zero: V(r 0 ) = 0 
a) Reference point: r 0 = R 0 (sea level) 


V(r) = mg R 2 0 (j- - = mgR 2 0 (~ - \ ) 

fto t Hq /to + h 
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where h = height above ground. 


= m 9 flg h 

{) R 0 (R 0 + h) 

If h << Rq V'(A) = mg h 
b) Reference Point: r 0 = oc 


r(r) = 


mg fig 
r 


P. S. 

The potential energy per unit mass is called the potential . 

Yii) = jR o 

m r 


V'(r) 


Example 2: Linear Spring Force 



F 


F = kx 


v = J a k x d x = X -kx 7 \ B A = ^ k {x 7 B - x\) 


If the spring is initially unstretched, the potential energy of the spring for the 
elongation x is: 
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2.4 Applications of D’Alembert’s Principle 


Example 1 

Dynamics of an accelerometer inside a rocket. 



M = rocket mass 
k = spring constant 
m = mass of acceleration 
r = displacement along sensitive axis 
F r = thrust force 
g = gravity 

= aerodynamic force 


Reference Frame is fixed in rocket with origin R 0 . 
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m a = —kr + mg - m[R 0 + w x (w x R) + (w x R) + 2(w x v)] 


where R = distance of m from origin 0. 
Dynamics of Rocket: 


M R 0 = A/g 4- F a + F r 


m a = —kr + mg — m 


g + 77 + 77 + wx(uJxR) + (u>xR) + 2(u> x v) 
MM . 


The gravity force mg cancels! 

NOTE : 

The (steady-state) surface of a liquid inside a rocket is perpendicular to the com- 
bined thrust and aerodynamic forces. 

Example 2 : 

Particle on turntable which rotates with uniform angular velocity u> (no friction). 


a = — 2 {uj x v) - u x (w x r) 
u) = wk r = xi + y] 


(1) x = 2 ojy + u 2 x 

(2) y = -2 uix -f u 2 y 


Introduce : 

z — x + i y complex variable 
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Multiply Eq. (2) by i and add to (1): 

z + 2 i u) 2 — w 2 z = 0 

Assume : z = z 0 e at 

Characteristic Equation: s 2 + 2 i uj s — u; 2 = (s + i uj) 2 = 0 
Si = —i uj S 2 = it 1 Repeated Root! 

z = (z Q + Z\ t)e 

Remember : Variation of Parameters. 

Initial Conditions: z(0) = 0 i(0) = v 0 (COMPLEX) 

z = v 0 t e~ iwt 


The particle moves radially outward with uniform velocity which is superimposed 
by clockwise angular velocity. 

The path of the particle in (x, y) plane is: 

Set = Xq — ► z = x 0 t(co s ut — i sin uii), 

X — Xq t COSCjf, 

y = —xqI sin ui 
Polar Coordinates : 


x = r cos(f> 
y — r sin <f> 
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x 2 + y 2 = *0 t 2 = r 2 — » r = i 0 t 


Let <f) = — u>t 


Kinetic Energy : 


r = 

Archimedian Spiral 


T = ^ m v 2 = ^ m(i 2 + j/ 2 ) — -m [(i 0 coswt — i 0 * w sin uit ) 2 
+(io sin + io* w cos wt) 2 | 


T = l m x 2 0 + l mx 2 t 2 u 2 


From Archimedian Spiral: t 2 = t? 

2“n 


1 1 


T = - m in + - m u> r 


2 _2 


The second term is equal to the work done by the centrifugal force. 


W = Fdr = f 

yo -'0 


Fdr = I m u 2 rdr = - muj 2 r 2 


The centrifugal force for u> = const is a conservative force. For u> = w(<), it is an 
irrotational force. 
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NO TE: 


The Coriolis force does no work because it is perpendicular to the velocity. It is a 
wattless force like the Lorentz force in electrodynamics [F = F v = — 2m(w x v) • v] 

P.S. 

If the centrifugal force had been neglected as a higher order term of u> then the path 
would be a circle (dashed line) in figure. 
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Historical Note : 

During the British-German naval battle of the Falkland Islands (about 50° latitude), 
the British gun shots landed almost one hundred yards to the left of the German 
ships, because the firing tables had been calculated for Britain’s northern latitude. 


NOTE : 

For a projectile on the surface of the earth, the effective angular velocity is u = 
w 0 sin <j> where <j> = latitude. 
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Example 3 : Buys-Ballot Lane (1817-1890) 



In a cyclone, the winds rotate about a center of 
low atmospheric pressure clockwise in the southern 
hemisphere and counterclockwise in the 
northern. 


Example 4 : Foucault Pendulum (1851) 

X 

x = 2 u> sin (f>y + (w 2 sin 2 <j>)x - g — 

y 

y = — 2u;sin <j)x + (w 2 sin 2 (f>)y — g - 


Introduce: z = x + iy and repeating the steps of example 2 yields: 

z + 2i u i + (j — u 2 )z = 0 


where u = u> sin <f> (<f> = latitude) 

Assume: z = Ae*‘ 

s 2 -t- 2 i u s + j = 0 


s i 



s 2 = (-u - 
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General Solution: 


a = A\t + A 2 e' 3 ‘ 
z = (Aifyft* + .4 2 e _, Vf') e - , “ t 



The complex vector z rotates with clockwise 
rotation in the northern hemisphere while 
performing a pendulous motion. The rotational 
angular velocity is wsin© and the pendulous 
motion has frequency of ui p = 


Example 5 : “Incense Swing” 

Pilgrims to Santiago de Camposlella, Spain, visit the shrine of St. James to burn 
incense. The incense and charcoal are held in a large silver brazier hung from 
the ceiling. The brazier is set swinging with a small amplitude, and then it is 
pumped by about six men until it is swinging through 180°. The swinging makes 
the charcoal burn energetically for the pilgrims. The pumping is the interesting 
part: they do it by shortening the rope by about a meter each time is passes through 
the vertical; they release the same amount of rope when the container reaches its 
maximum height. How does this shortening and lengthening of the rope increase 
the amplitude? 
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Each time the pendulum is suddenly shortened by 8 and then lengthened by the 
same amount at the extreme position. 

The motion of the pendulum is 


(j> + Wq sin <j> = 0 


where = g jl 


and 


cfy 2 = 2u>q(cos < j > — 


cos^o) 


Tension : T = m(<^> 2 / + g cos<f>) 


<f> = () : W’i = 2<5m[2w 2 (l - cos <j> 0 )l + g) 


<p = (f> 0 : W’ 2 = 28mg cos <f>o 
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Energy Increase: AE = W’i - H’ 2 = 6 mg 6 (1 - cosd> 0 ) . 
Example 6 : “Twirling Ice Skater”. 



A whirling particle of mass m is pulled 
in by a string toward a fixed center at 0. 


ma = F — m(u x r) — mu x (w x r) — 2m(u x v) 
a) Angular Motion : ( <f> - direction) 


mur + 2 mur = 0 

du _ f T dr 


n . . f u du „ f T dr 

ur + 2 ur = 0 => / — + 2 / — = 
Ju/q U) J r a V 

ln( — ) + ln( — ) 2 = 0 
w 0 r 0 


or 


2 2 

mur = mu;or 0 (Angular Momentum is conserved) 

b) Energy 

Work done against centrifugal force 

W = / Erfr = — m f u) 7 rdr — — m(u;orn) 2 f — 

^ A 0 r 3 


IT = -m(w 0 r 2 ) 2 ( - i) where r 0 > r. 


Difference of kinetic energy: 


1 


AT = T - r 0 = -m(wV - u; 2 r 2 ) 


AT _ J /KM 22 \ 1 

Ai - r m ( — w 0 r 0 1 = -m 


VO 2 > 

' 0 * c 0 2 2 

— — " 

r * j 


AT=im(w„r 0 ^(i-i)=ir 


Q.E.D. 
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Chapter 3 

Dynamics of a System of Particles 


The laws of motion will now be extended to a group of N particles which act on each 
other within a certain bounding envelope. The interaction forces between these 
mass particles are assumed to obey Newton’s third law: 


a;/, = -a;. 


3.1 Translation and Rotation 


Summing up all the forces (external and internal) Newton’s second law is: 


A r 

m t R, = F, + E R* fc 


*=i 


(3.1) 


where the sum extends over all N particles. Then summing over all N particles 
gives: 


E ”, a, = E f, + •£ £ R; t (3,2) 

Because of Newton’s third law the total interaction between the particles becomes 
zero. Therefore: 

E R-. = E F * = F ( 3 - 3 ) 
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where F is the total external force. 


Remember : 

R, is the position vector of the 

i-th particle in the fixed (inertial) frame. 

The associated acceleration R, is 
called the inertial acceleration. 

(R, is the inertial velocity.) 

Using Equation 1.4 of Section 1.4, Equation 3.3 can be written in the form: 

Y m, [Ro + (wxr,)+wx(«x r.) + 2 (w x v.) + a,] = F (3.4) 


This is the translation equation of motion referred to a moving frame. 

Very often (but not always) it is of advantage to select the origin of the moving 
frame such that: 


Y ™, r , = 0 


at all times 


(3.5) 


This means also that: 


Y v, = o 


and 


Y rn i a i = o 


(3.6) 


In other words, the origin of the moving frame is made to coincide with the center 
of mass (C.M.) location of the particle system. By definition the C.M. is: 


R, = — Y m i R-j where m = Y' 

m 


m, 


(3.7) 


In this case the translation equation (3.4) simplifies to: 


m R 0 = m K c = F 
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(3.8) 
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This is the center of mass law which states that the C.M. moves as if the sum of all 
external forces were acting on the total system mass concentrated at the C.M. 

NOTE : 

The total force F has to be determined for the actual system of particles and not 
for the total mass at the C.M. 

Examples : 

Exploding bomb shell; chair being thrown out of window. 

Next consider taking moments of Equation 3.1 about the origin 0 of a moving 
frame: 


r. x (m,R t ) = r, x (F, + £ H-ii.) 


(3.9) 


Again summing over all N particles yields: 


£ r, x (m,R,) = £ [r, x (F, + £ R'j] = £ r, ■ x F, + £ £ r, , x R'* (3.10) 


Because of Newton’s third law, it can be shown that the interaction effect again 
vanishes. 

Therefore: 


53 r, x (m, R t ) = 53 r > x = L o (311) 


where L 0 is the total external moment of the forces (torque) taken about the original 
O of the moving frame. 
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NOTE : 


The moment of a vector depends on the reference point of the position vector, i.e., 
changing the reference point changes the moment. 

equation 3.1 1 can also be expressed in terms of a moving reference frame via Equa- 
tion 1.4 of Section 1.4. The result is: 


Y, m . r, x [Ro + u x r, + u x (u> x r.) + 2 (u x v,) + a,] = L 


(3.12) 


This is the rotation equation of motion. 

If the origin of the moving frame is again made to coincide with the C.M. Equa- 
tion 3.12 becomes: 


E m i T x [w x r, • + u> x (w x r.) + 2(w x v.) + a,] = L 0 (3.13) 


Comparing Equation 3.8 and Equation 3.13 it is seen that choosing the C.M. as 
reference point O leads to a dynamic decoupling of the translational and rotational 
motion. 

The equations are however, in general, coupled through external forces like aerody- 
namic forces etc. The latter depends on the translational velocity and orientation 
(e.g., angle of attack) of the system. The forces (and therefore, the associated 
torques) are, in general, given as 

F = F(R, v, a,u>) (3.14) 


where symbolically, R — position, v = velocity, a = orientation, and u) = angular 
velocity of the system. 

The following special cases can be encountered: 

CASE A: F = F(R, v) 
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Solve : 


first translation 3.8 


then rotation 3.13 


CASE B: F = F (a, u) 

Solve : first rotation 3.13 then rotation 3.8 or 3.4 

3.2 Linear and Angular Momentum 


Revisiting the translational Equation 3.3 it is possible to reformulate it by introducing 
the concept of the linear momentum of a system. It is defined as: 


P = X m i = Y. “i v i 


(3.15) 


With it, Equation 3.3 becomes: 

P = F 


ra,- = CONSTANT (3.16) 


The time rate of change of the linear momentum of a system is equal to the total 
external force acting on it. In particular, if the external face is zero, the linear 
momentum of a system remains constant, i.e. 

P = P 0 for F = 0 (3.17) 


This is the conservation law of the linear momentum . 

A similar relationship can be established for the rotational equation by introducing 
the concept of the angular momentum (moment of linear momentum) of a system. 
It is defined as: 


Ho = Y. r « x ( m > H.) = 5Z r, x (™. v.) 


(3.18) 
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Taking the time derivative we obtain 



H 0 = Y, *«' x ( m < 

v,-) + 5Z r,(m,R,) 

(3.19) 

Because v t 

= v 0 + r,, it follows then that 


EJ5L 

From Eq. (3.11) — ♦ L 0 = Y 

r,- x (m,- R,). 



H = Y. m i r. 

x (v 0 + i\) + L 0 

(3.20) 

and finally: 





Ho = £ (m, 

r t ) x v 0 + L 0 

(3.21) 


The rate of change of the angular momentum about the origin 0 of the moving 
frame is equal to the total torque about this origin plus a term depending on the 
velocity of the origin and the center of mass of the system. This is slightly different 
from our previous find for the translational equation. However, the rate of change 
of the angular momentum is equal to the torque for the following cases: 

CASE A : Origin is fixed v 0 = 0 

CASE B : Origin is C.M. Y m, r, = 0 

CASE C : VollE m, r, = m r c 

where r c is the velocity of the mass center relative to the moving frame as seen by 
an inertial observer. NOTE : 

Some authors define the angular momentum of a system about the origin of a 
moving frame as: 


H 0 = ]T r, x (m, r : ) 


(3.22) 


where r, is the relative velocity of a particle as viewed by an inertial observer. 
Taking similar steps as above leads to the final result: 
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(3.23) 


H 0 = L 0 — X] m * r > x Ro = Lo - r, x (mR 0 ) 


Regardless of how the angular momentum is defined, the time rate of change of the 
angular momentum of the system is equal to the external torque when referred to 
the C.M. or a fixed point 0. In particular, if the external torque about these two 
points is zero, the angular momentum of a system remains constant: 


H = H 0 if L 0 = 0 


(3.24) 


This is the conservation of angular momentum . 

NOTE : 

Linear and Angular Momentum are very useful concepts when properly understood 
and applied. However, force arguments based upon Equation 3.4and Equation 3.12 
provide a more direct insight into the dynamic behavior of a system. 

Example 1 : 

A person stands in a boat and walks in the boat a certain distance and then stops. 
How far will the boat move? 



m B V B + m p (V B + v p ) = 0 


T T 

\ B — ; 

m B + m p 


Set 



and 


v p = 


d[ 

dt 
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m p l 

x B = 

m B 4- m p 


The final result does not depend on the time history of v p . 

Explain motion using force arguments via Equation 3.4: £ m,(R 0 + a,) = 0. 

P.S. : 

The common way to solve this problem is to use the principle of motion of the mass 
center (see Equation 3.8) which states, that in the absence of external face, the 
mass center will remain at rest or in uniform motion: 


M X C — (2J m i a: t )bEFORE — (5Z m i x i) AFTER 

TR p 4 " tbr S — Tn p 4 - ttib(S 4 ~ 

Set : i 2 = 4 xb T l 

TRpX\ 4" m B S = irip^Xi 4- Xb 4~ /) 4~ trb{S 4~ xb ) 
Solve for xb 


Example 2 : 

A horizontal circular disk can rotate freely about its vertical axis. Along a circular 
track of radius 1, a particle Q starts travelling with a constant speed v. What 
angular velocity will the disk acquire if it was initially at rest? What happens if 
the particle stops? 


H 0 = £ m, r, x v, = 0 

The velocity of the mass elements of 

the disk is v, = u x r, 

The angular momentum of the disk is then: 
Hi - (£ m, r 7 t )u) = / u> 
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where I = moment of inertia without small mass m. 


The angular momentum of Q is: 

h = m l(v + ul ) 


Hi + h = Iu> + m /( v + uil ) = 0 
ml v 

u j = 

I + mP 


When the particle stops, the disk stops too, but it has rotated through a finite 
angle <j>. 


3.3 Kinetic Energy and Work 

The concept of kinetic energy and work will now be extended to systems of particles. 
The kinetic energy is defined as: 

T = \ E m > ft? = \ E m .(k 0 + r ,) 2 
= \ E m * \ E m « r 2 + E m « r, • Ro (3.25) 

or 

T = - m R 2 + - YL m t r, 2 + m Ro • r c (3.26) 


Theorem of Konig 


For the case where R 0 = Rc, Equation 3.26 reduces to 

T =\ m &c + \ E m ^ 2 


(3.27) 
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To derive the expression for the total work done by all the forces acting on all the 
particles m, of the system it is assumed that R 0 — Rc« 


H - £ H’i = E [*' ( F - + Z R W ( dR c + dr,) (3.28) 

^ a i 

or 

H- = E f BC (F, + ER'j) ■ dR c + Y. f < F ' + E R i) ' (3.29) 

■/.4c ^ -4 , 

Because of Newton's Law this simplifies to: 

w = [ Bc F . rfR c + X f B ' (F, + X R*,) • dr, (3.30) 

J Ac J Ay 

Now for each particle, the principle of work and kinetic energy applies. Therefore, 
using Newton's third law, we get 

= \ m R c \a C c + \ E m ' r ? (3.31) 

Comparing Equations 3.28 and 3.30 it follows from the center of mass law (Equa- 
tion 3.8) that 

j Bc F-fficTmR’ |£ (3.32) 

J Ac *■ 

and 

E J°' (Pi + R o) ' = E r? If; (3.33) 

Therefore the work done by the external and internal forces is equal to the increase 
in the kinetic energy of relative motion. The velocities r t relative to the C.M. can 
arise from rigid body rotations in which the distances between the particles do not 
change, as well as the more obvious case of changing particle separations. 

In cases where the external and also the internal forces can be derived from a 
potential energy, the total energy of the system is again conserved. 

E = T a + \ A = Tb + V b (3.34) 
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3.4 Variable Mass 


Occasionally dynamic systems have to be dealt with whose mass varies with time. This 
is caused by mass particles either leaving or entering a certain boundary envelope 
(control volume). Consider for the following preliminary example. 


V 



Two bullets are fired from a vehicle which is rolling on rails without friction. Cal- 
culate the velocity of the vehicle after the bullets have been fired (a) simultaneously 
or (b) in sequence. The muzzle velocity is C relative to the barrel and M is the 
vehicle mass without the bullets. Applying the principle of linear momentum we 
obtain 


(a) M \\ + 2 m(V'i - C) = 0 l\ = 


2mC 
M + 2m 


After first bullet is fired 

(b) 


V? = 


mC 


M + 2m 


After second bullet is fired 


1-2 = m C(T7” h 77 - 0 — ) 

v M + m M + 2m 


Note that l j > V). Consider now the case where N bullets of mass m are fired 
sequentially: 


V = C 


£ 


A=1 


m 

M + km 
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If the bullet size becomes infinitesimal small and N —> oo the summation can be 
replaced by integration to yield: 

rMo dm 

V = c I — 

Jm m 

where Mq is the total initial mass including the bullets, 
and finally: 


V = C In 


Mo 

M 


This is the well-known velocity equation for a rocket in free space (no gravitational 
and aerodynamic forces). The velocity is seen to be proportional to the muzzle 
velocity (exhaust velocity). An interesting step can be taken by differentiating the 
velocity equation with respect to time: 


or 




M f = —C M = F t 


Ft — thrust force. 

This is the basic force equation for a rocket and thrust force 

F t = -C M 

Notice that the fuel flow M is negative because the system is losing mass. 

A more general approach to the dynamics of a variable mass system can be taken 
by examining the translational Equation 3.4 and the rotational Equation 3.12 

A) Translation (Thrust force) 

To gain an understanding of the thrust producing mechanism we set u = u = 0 
and assume also that there are no external forces acting on the system. (F = 0). 
Equation 3.4 reduces then to: 

m R 0 + ]T m, a , = 0 m - m 0 + J2 (3.35) 
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where m is the total mass of the system at any point in time and m 0 represents 
the mass of the main body which remains constant. The origin 0 of the moving 
reference frame is assumed to be fixed in the main body. It is obvious that the 
thrust producing mechanism must originate from the second term of Equation 3.35. 
To show this, we introduce a modification in the notation to rewrite this term as: 


X ] m i a > 


lim 

A(—0 


Z Am.Av, 
A i 


(3.36) 


Consider now the case in which an abrupt change in velocity Av is imparted to a 
large number of very small particles. Then Equation 3.35 becomes 


dmi 


£ m, a, = Y. A v . = T. Av . 


(3.37) 


Notice that m, is the mass per unit time which undergoes the velocity change Av, 
relative to the moving frame. For steady state conditions this becomes the rate at 
which mass leaves the system. The acceleration which the main body attains is 
occurring during the short acceleration phase of the moving mass particles. It is 
important to realize that this acceleration is not caused by the fact that mass is 
leaving or entering the system. From Equation 3.35 we can now write: 

m R 0 = - 5 Z "i. Av, = Ft (thrust force) (3.38) 

It is seen that the term (-mAv) acts as an effective force on the main body. By 
the nature of the above derivation the term m is always to be taken as positive 
regardless of whether mass is leaving or entering the system. The difference be- 
tween the two is that particles leaving the system are experiencing a rapid relative 
acceleration, whereas particles entering the system are rapidly decelerated. Both of 
these processes result in a finite change Av of the relative velocity. 

In a rocket the thrust is produced by accelerating gas particles from zero relative 
velocity rearward until they reach a relative exhaust velocity C. The magnitude of 
the thrust is therefore, 

F t = rhC (3.39) 


It is common rocket engine practice in the US to characterize the performance of 
an engine by the specific impulse (specific thrust) which is defined as: 

/, - — = -El— — — seconds (3.40) 

W 17100 9o 
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The specific impulse is the thrust per propellant flow rate. It is really a measure 
of the effective exhaust velocity C and its only merit is that it has the same unit 
in the metric as in the customary system of units. Typical values of the specific 
impulse are in the range of 200-450 seconds. 

Space Shuttle SSME : I s = 410 seconds (average) 

Space Shuttle SRB : 1$ = 265 seconds 

The effective exhaust velocity is simply C — I s g 0 ■ In the metric system g 0 = 9.81 
m/s 2 ~ 10 m/s 2 and the effective exhaust velocity is just ten times the specific 
impulse. 

P.S. 

The speed of sound in air is approximately 300 m/s. It is often useful to express 
velocities in terms of the Mach number. As an example, the specific impulse Is 
= 410 seconds corresponds to an exhaust velocity of approximately 14 Mach. The 
speed of light is C = 300,000 km/s = 3 x 10 8 m/s. It is equal to 10 6 Mach (1 
Megamach!) 

A very important case of mass flow is that of steady fluid flow in pipes. The 
continuity condition requires that the flow rate m be constant, i.e. 

m = p AiVi = pA 2 v 2 (kg/s) (3.41) 


where Ax is the entrance cross section and A 2 the exit cross section, and Vi and v 2 
the corresponding velocities. 
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According to Equation 3.38 the force exerted by the liquid flow on the pipe is: 

Ffl = -m(v 2 - V)) Euler’s equation (3.42) 


B) Rotation 

If the dynamic system has an angular velocity internal mass flow generation pro- 
duces an additional effect coming from the Coriolis term of Equation 3.4. Concen- 
trating only on this effect, the acceleration of the main body is: 

m R 0 + 2 51 m ii u * v i) = 0 (3.43) 

The Coriolis term can be modified again as previously to: 

, \ „ V"' A m i( W X ^ r ' ) (o 1 A 

2£m i (wxvj) = 2 2 </ {tH) 


Taking the limits A t — ♦ di and Am, 
as: 


2 5Z ™.( w x v >) = 


— > dm,, Equation 3.43 can formally be written 
— x Ar,) = 2 x Ar ') (3-45) 


The physical interpretation is that it represents the effect of an internal flow rate 
which extends over a distance Ar. The resultant acceleration of the main body is 


therefore: 

m R 0 = -2 ^ m,(u> x Ar,) (3.46) 


In general, the flow rate will extend over a finite length, such that the total effective 
force of the flow rate will be obtained by integration. Therefore, we obtain finally: 

m R 0 = -2 u) x 5Z / ™«^ r > (3-47) 

J r A 


where the integration extends over the stream lines. The effective Coriolis force 
stemming from internal flow rates is seen to be dependent on the geometry of the 
flow. It is therefore, quite difficult to make a statement as to its overall effect on 
the system acceleration. 
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These two effective forces caused by mass flow rates affect, of course, also the 
rotational motion as governed by Equation 3.12 by producing concomitant effective 
torques. Of particular interest is the effect of the Coriolis torque caused by a flow 
rate m because it gives rise to the so called jet damping effect. A simplified example 
is given to illustrate the situation. A rocket rotates about a transverse axis through 
the center of mass at 0. 



A uniform mass flow rate is assumed to exist along the x-axis only extending from 
x — l\ to x = / 2 which is the nozzle exit. 


The Coriolis torque is: 

Lc = — 2 J m r x (u> x dr) 
and with u) = wk and r = x i: 


L c — —2 mu — - — — = -mw(/j — /^) 


(3.48) 

(3.49) 


A few special cases are worthy of note: 

a) l\ = ^ Lf = 0 

b) l\ = 0 L c = -mwlj 

c) / 1 = — l? L c r 0 
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Case (b) is often given in textbooks as the general term for the jet damping using 
erroneous angular momentum arguments. To understand the jet damping effect, 
recall the example of the twirling ice skater. 

NOTE : 

The rotational equation of motion is 

Id) = — /j) 


It is seen that the Coriolis torque is causing an angular deceleration proportional 
to the angular velocity u> of the rocket. This can be physically interpreted as an 
effective damping . Since the Coriolis force is wattless, the energy dissipation has to 
be caused by the centrifugal force. 


3.5 Impact Dynamics 


During the impact of two bodies, very large forces act for a very short time. Such 
forces are called impulsive . Because they are so large, all other forces (e.g. gravity) 
can be neglected in their presence. When impulsive forces act on a body, the 
velocities undergo an instantaneous finite change (Av^O) whereas its position and 
orientation remain unchanged (Ar = 0). The impulsive forces are specified by their 
short duration time integral: 


F = 



F di 


6 = short time 


(3.50) 


They are treated similarly to the Dirac delta function. The linear and angular 
velocity changes of a body during impact are obtained by integrating the equations 
of motion 3.4 and 3.12 with respect to time. Since only rigid bodies will be 
considered here, the relative velocity v and the relative acceleration a are zero. 
The resultant impulse equations are algebraic equations with the velocity changes 
as unknowns. 


mAVo + Aw x 


(mr c ) = j 


Fdt = F 


(3.51) 
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(mr c ) x AV 0 + £ ™.r t x (Au> x r t ) = / L 0 <ft = L 0 


(3.52) 


These are two algebraic equations which have to be solved for the velocity changes 
AV 0 and Au>. However, another equation is required because the impulsive force F 
is also unknown. (The impulsive torque L 0 can be determined from the impulsive F 
by taking its moment about the origin O). There will, of course, be a set of impact 
equations for each of the two colliding bodies. 

The additional equation sets up a relationship between the normal components of 
the relative velocities of approach and separation of the two bodies. 


V 2M ~ Ul N 
e = 

u lN — ^2/V 


(0 <e < 1 ) 


(3.53) 


where e is called the coefficient of restitution . It is notationally convenient to use 
the letter u for the velocity before impact and v for the velocity after impact. The 
coefficient of restitution depends on the material of the colliding bodies, on their 
geometries and also upon the impact velocity. All impact conditions lie between 
the two extremes e = 0 and e = 1. 

Case A : e = 1 (elastic) 

According to Equation 3.53 the velocity of approach is equal to the velocity of 
separation. No energy loss. 

Case B : e = 0 (plastic) 

The two bodies stick together after impact. Maximum energy loss. 

Impact phenomena are almost always accompanied by energy losses. The higher the 
impact velocity the more energy loss occurs. The coefficient of restitution actually 
approaches UNITY (no energy loss) as the impact velocity goes to zero. Energy 
is lost through heat generation, generation and dissipation of internal vibrations 
(elastic stress waves) and sound energy. 

NOTE : 

The impact dynamics equations are a set of linear algebraic equations. 
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Example 1 : 




Vl 2 



Consider first the central collision of two spheres. 


F = 


— - ui) = m 2 ( v 2 - u 2 ) 


1) Plastic Collision : 


rriiUi + m 2 tt 2 = + m 2 i; 2 


v x = t» 2 = v 


rriiUi + m 2 u 2 

v = 

+ m 2 


2) Elastic Collision : 


mxUi -I- m 2 u 2 = mill! -f m 2 u 2 MOMENTUM 


m!«j 4- m 2 i/2 = myv"l + m 2 U2 ENERGY 


mi(v! - ui) = — m 2 (n 2 - u 2 ) 
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TTli(v J — Uj) = — m 2 ( — U 2 ) 


i?l + If] — 1^2 ”1" ^2 


or 


(v 2 - Ui) = («1 - u 2 ) 

SEPARATION VELOCITY = APPROACH VELOCITY: 1 = 

«1 -»2 


Solve for the velocities after impact: 


v x =Ui- 


2 m 2 

m i + m 2 


- u 2 ) 


2mj , 

ii 2 — u 2 + ( u 1 

mi + m 2 


- u 2 ) 


3) Inelastic Collision : 


v 2 - «i _ SEPARATION VELOCITY 
u l -u 2 ~ APPROACH VELOCITY 


miiii + m 2 i >2 = rriiUi + m 2 u 2 


Vi - V 2 = —e(ui - a 2 ) 




= «i - 


m 2 (l + e) 
m l + m 2 


(«i 


— u 2 ) 
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v 2 = U 2 + 


m^l + e) 
mi + m 2 


(ui - u 2 ) 


Loss of kinetic energy : 

A T = ^ [m^u 2 - uf) + m 2 {u\ ~ ^ 2 )] 


After some algebraic manipulation the final result is: 



m 1 m 2 

mi + m 2 


(1 - e 2 )(ui - u 2 ) 2 


for e = 0 : AT = — u 2 ) 2 


where - = — + u = reduced mass, a = Greek “mu” 

m 1 m2 r 

Example 2: Ballistic Pendulum 



For m 1 : F = — u() 

For m 2 : F = m 2 (t4 — t> 2 ) 


Combining Equations 3.54 and 3.55 yields the conservation law of linear 
turn. 

mi Ui + m 2 u 2 = mi v[ + m 2 t>2 


(3.54) 

(3.55) 


momen- 

(3.56) 
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The coefficient of restitution relation is: 


v' 2 -v[ SEPARATION VELOCITY 
v 1 -v a ~ APPROACH VELOCITY 


Initial condition for m 2 : u 2 = 0 


We obtain for the unknowns v[ and i> 2 : 


. mi — m 2 e 

Vi = «i 

mi + m 2 


/ (l + e)mi 
'2 = ; wi 

mi + m 2 


The angular velocity after impact is 


Example 3 


a> _ /; _ (! + e ) m i«i 

0 — Vo/l = TT 

(m, + m 2 )/ 


Mass m, moving along the x-axis with velocity v, hits m 2 and sticks to il 
three particles are of equal mass m, and if m 2 and m 3 are connected by 
massless rod, find the motion after impact. 



(3.57) 


. If all 
a rigid 


1 



Linear momentum in normal direction: 


m i>o 


= m v' + 2 m x 


Linear momentum in tangential direction: 

0 = m v' y + 2 m y' 


Angular momentum about C.M. of dumbbell: 


. / , U 0 ]/ / ^ \ V 'x fl\ 

^ ^ + - m (o) 


‘2 y /2 


JjL 

V2 


Coefficient of restitution: 


(&+*•)-< 

e — 

Case A : Sliding Rebound (/i = coefficient of friction) 

R = fi F = -m v' y 

F = m{y o — v x ) FOR PARTICLE 

Combining these two equations yields: 

/*(»0 - * 4 ) = ~ v y 


Rearranging the last five equations yields 

(P) v' x + 2i' = v 0 

(2') v' y + 2y = Q 

(3') v' x + i’t n/2 - i>; = u 0 


(3.58) 


(3.59) 


(3.60) 


(3.61) 


(3.62) 
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~ V 'r + + *' = e ^0 


(5') n v' x - v' - n v 0 


Solving these equations for the unknowns results in: 


, _ (3 -fi - 4e) 

(7 - n) V ' 


, -4/i(l + e) 


, 2(1 + e) 

1= 


•/ _ Ml +£) 

* (7 - #.) "° 


_ 2y/2(l - p)(l + e)t> 0 

(7 - /0* 


Case B : Normal Rebound 

Equation 3.62 is replaced by normal rebound condition: 

(5") < = v'-£k ' 

All other equations remain the same. Solving for the unknowns: 


5 - 7e 


1 + e 
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i' = ^ (1 + e) vo 

y‘ = h C 1 + e ) w o 



(1 + e) vp 

2s/l 


The coefficient of friction above which normal rebound occurs is obtained from 
calculating the two impulsive forces: 

R = fi f = -m v' y = m u 0 

7 

F = m(v o - v' x ) = —(1 +e)mv 0 
^ = 2/7 


Example 4 : 

A sphere m\ is projected horizontally against a carriage m 2 w’hich is backed up by 
a spring k. If the coefficient of restitution is e and the surface is perfectly smooth 
(fi = 0). Calculate the rebound velocity v[, the rebound angle 6 and the maximum 
travel 8 of the carriage after impact. 



The prime indicates the condition after impact. 
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Impulse Equations : 


m i (x[ — ij) = —F sin a 

(3.63) 

m i Oil ~Vi) = F cos a 

(3.64) 

m 2 (i 2 — i 2 ) = F sin a 

(3.65) 

y 2 = 0 carriage is on rail 

(3.66) 


Normal Velocities: 


v in == i\ sin a 

(3.67) 

l; lN = sin Q — cos Q 

(3.68) 

^2A f = 0 

(3.69) 

t’jjy = sin a 

(3.70) 


Restitution Equation : 

_ l: 2 \ ~ v[iv _ z 2 s ^ n Q ~ (^1 s » n ° ~ y[ sin a) 
Via- — x sin a — 0 


(3.71) 


Unknowns : 

* 1 > Vli ^ 2 : F 1 V?' 

The following steps are taken: 

A) Eliminate F : 

1) Combine Equations 3.63 and 3.65: 

x[ + m 2 x j = m ] x l linear momentum for m l in x-direction (3.72) 
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2) Multiply Equation 3.63 by cos a and Equation 3.64 by sin a: 


mi cos a + mi y[ sin a = miii cos a 


(3.73) 


The linear momentum along the smooth surface is conserved for m^ 
B) Solve Equations 3.71, 3.72, 3.73 for x[, y[, and x' 2 


C) 

Determine fl: 

• / 




* a Vl 
tan v — — 

x[ 

(3.71) 

D) 

Determine i/: 

v' = J(i\y + (!i !) J 

(3.75) 

E) 

Determine fc \ (x^f 

-i 



Example 5 : 

A uniform bar of length £ and mass m falls on a horizontal floor with velocity 
v 0 = u 0 i — uoj. The bar falls without rotation. If a coefficient of restitution e 
and a coefficient of friction /i exists between the floor and the bar, determine the 
minimum friction coefficient for a normal rebound and the velocities after a sliding 
rebound. 

y 



Impulse Equations: 


F = m(y + v q) 


(3.76) 
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R = —m (x — u 0 ) 

F- cos 8 = lui 
■ £ 

R - sin 8 = /w 2 


The total angular velocity of the bar after impact is 


UJ — UJi — U) 2 


Restitution Equation : 


e — 


_ y + | cos 8 (wi — >jj 2 ) 

Vo 


A) NORMAL REBOUND : (R < fiF ) 

£ 

i + (- sin 8) (u>i - w 2 ) = 0 

Substituting Equation 3.76 in 3.78 and 3.77 in 3.79: 

, . , t m£ 2 

m {y + ^o) - cos ^ = — ■ w i 


or 


Wi = 


2 12 

6 (y T n 0 ) cos 61 

7 


/ • a 

m{Uo — x j - sin 0 = u> 2 


or 


— 


6(uq — i) sin 8 


(3.77) 

(3.78) 

(3.79) 


(3.80) 


(3.81) 

(3.82) 

(3.83) 

(3.84) 


Inserting 3.83 and 3.84 into Equation 3.80 and 3.81: 

e w 0 = y + 3 cos 8 [(y + v 0 ) cos 8 — (u 0 — x) sin 8] (3.85) 
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x = 3 sin 8 [(u 0 - x) sin 8 — (y + v 0 ) cos 8] 


(3.86) 


(3 sin 6 cos 9)x + (1 4- 3 cos 2 9) y = e v o 
— 3 v 0 cos 2 9 + 3u 0 sin 9 cos 9 


(1+3 sin 2 8) x + (3 sin 9 cos 9) y = 3 sin 2 8u 0 - 3v 0 sin 8 cos 9 
Let 9 = 45° : 

cos 8 = sin 8 = = 

V2 

3 . 5 . , 3. 3 

2 x + 2 y = ( e - 2^ V ° + 2 “° 

5.3. .3. 3 

- x + - y = u 0 (-) - - v 0 


Final results for velocities after impact are: 


(5e — 3) v 0 + 3 uq . 3 {uo — (1 + e)v 0 } 

y= § * = 8 


(jJ = UJ\ — UJ2 — 

The impulsive forces are: 

R = m [ 


{2(1 + e) v 0 - u 0 } 


F — m 


2^21 

\ 5u 0 + (1 + e) Vo 1 

l 8 J 

3u 0 + 5(1 + e) Vo 1 

8 1 


Normal Re bound Condition: R < fiF 


(3.87) 

(3.88) 

(3.89) 

(3.90) 

(3.91) 

(3.92) 

(3.93) 

(3.94) 

(3.95) 
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3 /i 0 + 5(1 4- e) v 0 
5 //q + (1 + e) w 0 


SLIDING REBOUND : R = n F 


R = fi m(y + v 0 ) from Equation 3.76 


From Equations 3.77, 3.78, and 3.79: 


H rn (y + v 0 ) — —m(x — u 0 ) = rn(u 0 — x) 


From Eq. (5): 


i ■ \ ^ mP 

H rn (y + v 0 ) - sin 9 = -yy- 

/ • v. t- a mP 
m [y + v 0 )~ cos 9 = -yy Ul] 

6(?/ + v 0 ) cos 9 6 fi(y + u 0 ) sin 9 

u>i — u > 2 = - 

l £ 


evo = y + 2> cos 0 [(y + u 0 ) cos 9 - /i (y + v 0 ) sin 9) 


e v o — y [1 + 3 cos 0(cos 9 — fi sin 0)] + 3u 0 cos 6 (cos 9 — fi sin 9) 


Let 9 - 45° : 


Vo_ e — 3 cos 9 (cos 9 — fi sin 0)] 
1 + 3 cos 9 (cos 9 — n sin 9) 


= {2e - 3(1 - y)} . _ 2 /i(l + e)v 0 

2 + 3(1 — (i) X ~ Uo 2 + 3(1 — p) 

* = ^-* 2 = o 

v^[2 + 3(l 


68 


(3.96) 

(3.97) 

(3.98) 

(3.99) 

(3.100) 

(3.101) 

(3.102) 

(3.103) 

(3.104) 

(3.105) 

(3.106) 
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Chapter 4 

Dynamics of a Rigid Body 


This section is almost entirely treating the rotational dynamics of a rigid body. Later 
we will also include internal moving parts whose movements are often prescribed 
time functions. First the general equations of motion of Equation 3.12 are brought 
in a more suitable form. 

4.1 Euler’s Equations 

For a rigid body the relative velocity v and the relative acceleration a in Equa- 
tion 3.4 and Equation 3.12 of the preceding sections are zero. For the case of the 
general motion of a rigid body, it is customary to choose the center of man as the 
origin of the reference frame and have the reference frame fixed in the body (Body- 
fixed reference frame). If there is no external coupling between the translation and 
rotation, the stational motion can be analyzed separately. 

With these assumptions Equation 3.12 is: 

£ m, r, x (w x r.) + E m, r, x [w x (w x r,)] = L 0 (4.1) 

The first term is seen to be the (negative) Euler torque, whereas the second term 
is the (negative) centrifugal torque. This equation can be rewritten b^ using the 
following vector identity: 

rx[wx(wxr)]=wx[rx(ux r)] (-L2) 

The rotational Equation 4.1 becomes then: 

E m, r, x (tii x r() + u x E mi r, x (w x r.) = L 0 (4.3) 
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At this point it would be possible to introduce an orthogonal reference frame and 
convert Equation 4.3 to three scalar first order differential equations for dynamics 
analysis and computer programming. However, it is possible and very useful to em- 
ploy a notation known as vector-dyadic notation which, unlike the matrix notation, 
is independent of a reference frame in the same way as vectors are. 


4.2 Vector-Dyadic/Matrix Notation 


A new vector operation is introduced which is formed by the juxtaposition (side 
by side storage) of two vectors for the purpose of later taking scalar and vector 
products with an ordinary vector. This vector operator is called the dyadic product 
and is given by: ~ 


P = AB (4.4) 

where A is called the antecedent and B is the consequent . The dyadic product is 
not commutative because AB / BA. However, the distributive law holds: 

A(B, + B 2 ) = AB t + AB 2 

(Ax + A 2 )B = A y B + A 2 B (4.5) 

The juxtaposition of two vectors AB is also called a dyad. 

NOTE : 

Dyads will be designated by capital script letters. 

The sum of dyadic products is called a dyadic : 

V = AiB t + A 2 B 2 + . . . A„B„ (4.6) 

The dyadic obtained by interchanging the order of the A, and the B, is called the 
conjugate of V: 

V c = Bj Aj + B 2 A 2 + . . . B n A n (4.7) 

The dyadic which is equal to its conjugate is called self-conjugate or symmetric. 

Because the distributive law holds for a dyadic product any dyad AB can be written 
in form of a dyadic as follows. 
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Expressing the vectors A and B in terms of orthogonal unit vectors e 1 ,e 2 ,e 3 , we 
get: 


V = AB 


= (Ajei + A 2 e 2 + A3e3)(JBiei + B 2 e 2 + £3^3) 
= Aj^ieiei + Aif? 2 e 1 e 2 + AiB^^e^ 

4 - ^2^i e 2 e i + A 2 £ 2 e 2 e 2 + A2B3G2G3 
+ + A 3 J 9 2 e 3 e 2 + 


(4.8) 


This is called the nonion form of the dyad because it contains nine components. 


The following rules apply for the scalar and vector product of a vector with a 
dyad: 


V R = (AB) • R = A(B • R) 

Vector 


R • Z> = R (AB) = (R- A)B 

Vector 


V x R = (AB) x R = A(B x R) 

Dyadic 


RxD = Rx (AB) = (R x A)B 

Dyadic 

(4.9) 

The unit (identity) dyadic is defined as 

£ — e^i 4* ©2 e 2 + e 3 e 3 


(4.10) 


Indeed: A • £ = S • A = A 


We are now ready to go back to Equation 4.3 and cast it in vector-dyadic notation. 
Observing the above rule of forming the scalar product of a vector with a dyadic 
we obtain: 


£ m,T, x (u> x r t ) = £ m, [(r,- • r,)u> - (r, • u>)r,] 

= £ m, [r?£ - (r,r,)] • u = I • u) (4.11) 

where 

I = £ m,- [r, 2 Z - (r,r,)] 

I = inertia dyadic 

The vector-dyadic form of Equation 4.1 is therefore 

I'W + wxI'W = L 0 (4-12) 
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It should be reiterated that this form of Euler’s dynamical equations is independent 
of a reference frame. 

To convert Equation 4.12 into scalar form we introduce an orthogonal body-fixed 
reference frame (body axes) whose origin is fixed in inertial space or coincides with 
the mass center of the rigid body. The position of a mass particle m, is then given 
by: 


r,- = ar.-ei + y,e 2 -I- z,e 3 (4.13) 

and the angular velocity of the rigid body: 

w = u> x e i -I- o> y e 2 + w*e 3 (4-14) 

The ensuing scalar equations can be concisely expressed by using matrix notation . 

To this end, define the moments of inertia as: 

hx = £ 771 ,( 2 /? + *?) 

Iyy — £ 771,(1? + Z?) (4.15) 

I zz = £ 777, (i? + y?) 

Also define the products of inertia as: 

I xy Jyx ^ 

Jxz $ zx — ^ TTl{X{Z{ (4.16) 

Iy Z I Z y — £ rriiyiZi 

K-Q.TE; 

For continuous bodies the summation is replaced by appropriate integration e.g.: 


/ TT 


j (y 2 + z 2 )dm 


etc. 
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Caution: 


Some authors define the products of inertia of Equation 4.16 using the opposite sign 
such that: 

I xy = £ m.x.y, etc. 

The moments and products of inertia can be arranged in matrix form. 


1 = 


I xx f xy 1 xz 
lyx f yy I xz 
I zx I xy I zz 


I = inertia matrix 

Any vector product between two vectors can be cast in matrix form: 

c = a x b 

Introduce the skew-symmetric matrix 

0 fly 

CL 2 0 flj 

dy flj 6 

Then the matrix notation for the vector product is 

c = ah 

where b and c are now column matrices. 


a = 


a = singular 


( 4 . 17 ) 


(4.18) 


NOTE: 

No confusion is likely to occur by designating a column matrix by boldfacing the 
letter, which is the same convention used for vectors, because it should be obvious 
from the context whether an equation is written in matrix form or in vector-dyadic 
form. 

The matrix form of Euler’s dynamical equations is then 

Iu) + CjIu) = Lq (4-19) 

It should be mentioned that in the dynamics area w'here only orthogonal (rotational) 
transformations are used a dyadic is essential identical to a matrix. For every vector- 
dyadic equation there exists an equivalent matrix equation. The advantage of the 
vector-dyadic equation is, of course, that it is independent of the reference frame. 
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Examples: 

A) Angular Momentum of a Rigid Body 


H = £ r, x m,v, = £ r, x m,(v 0 + u> x r,) ( 4 . 20 ) 

If the origin O is fixed in inertial space or coincides with the mass center we have: 

H = £ m,r t x (w x r,) ( 4 . 21 ) 

or 

H = J • us 1 = inertia dyadic (4.22) 

The corresponding matrix equation is: 

H = /u7 (4.23) 


B) Kinetic Energy of a Rigid Body: 

T = 2 S m t( v . v «) 

= 5 £ m,(v 0 + w x r,) 2 

= j 777 u 2 + v 0 • (w x £ m, r,) + ^ E m,(u? x r ,) 2 


If the origin is fixed or at mass center, the middle term becomes zero. 

The total kinetic energy is then composed of translational and rotational kinetic 
energy. 


T = T t + T r 

The expression (w x r,) 2 can be transformed as: 


(u> x r,) 2 = (w x r) (w x r,) 

= ■ [r, x (a> x r,)] 

Therefore, the rotational kinetic energy of a rigid body is: 

T t = i u) • I • u> 
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(4.24) 


(4.25) 
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The corresponding matrix equation is: 


T t = i u> T I u> (4.26) 

where the superscript T denotes the transpose of a matrix. 

4.3 Orientation Kinematics 


An important aspect of rotational motion in space is the specification of the relative 
orientation between two reference frames (coordinate systems). In many applica- 
tions the orientation to be defined is that between a moving reference frame (or rigid 
body) and a space-fixed (inertial) frame. There exist practically three schemes for 
doing this: 

a) Direction Cosines 

b) Euler Angles 

c) Quaternoins (Euler parameters) 

Other methods are sometimes used, but they are of little or no advantage over these 
three. 

As a general rule the unit vectors of an inertial frame will be designated by i,j,k, 
and those of a moving frame by ei,e 2 , e 3 . 

a) Direction Cosines 

In this method each unit vector of the 
moving frame e, is defined 
in terms of the three angles each makes with 
the three axes of the inertial frame. They are 
obtained by the corresponding scalar products 
ei • i = cos Qn — A n 
ei • j = cosoi2 = A n 
ej • k = cosa 13 = Ai 3 etc. 

The coefficients .4,* are called the direction cosines . The first index relates to the 
unit vector of the moving frame and the second to that of the inertial frame. 


Z 
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There are nine direction cosines and they completely define the relative orientation 
of two reference frames. We can w'rite in vector form: 


e i — An i + A12 j + An k 


e 2 — A%\ i + A 22 j + A23 k 


e 3 — i + A32 j 4 - T33 k 

It can be easily seen that the nine direction cosines are not independent. The 
relationship between them can be obtained by observing the fact that the unit 
vectors e, are of unit length and mutually orthogonal. Therefore the following six 
equations hold: 


ei • e x = 1 e 2 ■ e 2 = 1 e 3 • e 3 = 1 

• e 2 = 0 e r ■ e 3 = 0 e 2 • e 3 = 0 

These yield six relationships between the .4,* which can be compactly written as: 

E A » = s* 

w'here 5,j is the Kronecker Delta. This equation is known as the orthonormality 
condition of the direction cosines. 

A more convenient form is obtained by using matrix notation. Introducing the 
direction cosine matrix : 


An An A13 

A — | A 2 1 A 22 A 2 3 

A 3 i A32 A33 


the above orthonormality condition is then 


A A t = E 

where E = unit matrix. 

Premultiplying Equation 4.27 by A~ l yields also: 

A' 1 = A t 
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(4.27) 


(4.28) 
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Matrices having the property that the inverse matrix is equal to the transposed 
matrix are called orthogonal (rotational) matrices . There are nine more (not inde- 
pendent) relationships among the elements of A given by: 


ei x e 2 = e 3 , e 3 x ei = e 2 , e 2 x e 3 — e x 


They can be summarized using the adjoint matrix adj A as 


A t = adj A 


NOTE : 

In many dynamics problems the moving reference frame originally coincides with 
' the inertial frame. It is, therefore, often customary to refer to the moving frame as 
the “ new ” frame and the inertial frame as the “old” frame. Using this terminology 
one would say that the direction cosines express the new unit vectors in terms of 
the old unit vectors. 

It is also necessary to know how the components of a vector change from one 
reference frame to the other. 


Consider the arbitrary vector 

R = a: e! -I- y e 2 + z e 3 = A" i + V j + Z k 


Introducing the unit vectors e!,e 2 and e 3 yields: 


R = (Ani + d 2 i!/ + ^3i 2 )i 
+ (Ai 2 x + A 2i y + A 32 z ) j 
+ (Ai 3 i + A 23 j/ 4- A 33 z)k = A i + > j + Z k 

Equating components of both sides and using matrix notation we obtain 


A 

Y 

Z 


= a t 


x 

y 

z 


or X = A t x 


( 4 . 29 ) 


This equation expresses the old components of a vector in terms of the new com- 
ponents. To reverse the situation we solve for the new components by using the 
orthonormality condition of the rotational matrix A and obtain 


77 


X ' 

Y or x = ,4X 
Z 



We see that the components of a vector transform exactly like their unit vectors. 

Since the orientation of the moving reference frame is constantly changing with 
time the direction cosine matrix is a function of time. It is possible to derive a first 
order matrix differential equation by determining the time rate of change of the 
unit vectors ei,e2 and e 3 relative to inertial space. 

Employing results from Section 1.3 we can write: 

e i =u x d = A n i + ii 2 j+ 4 13 k 

e 2 = u> x e 2 = A?i i + A?? j + A^k 

= u? x e 3 = -4 31 i + A^ j + ^4 33 k 

Expressing the angular velocity of the reference frame in terms of the moving ref- 
erence frame we obtain 


u) — ufi e 2 + u> 3 e 3 

The above vector products can then be written as follows: 


(4.30) 


u,xe i ~ (^i e i + e 2 + e 3 ) x e x = -a> 2 e 3 + e 2 (4.31) 

= -w 2 (-4 31 i + -4 32 j + /4 33 k) + w 3 (.4 21 i + .4 22 j + .4 23 k) (4.32) 

The other cross products can be transformed in like manner. Equating coefficients 
and solving for the time derivatives of the direction cosines we arrive at the following 
compact matrix differential equation: 


^ ^ A kinematical (Poisson) differential equations (4.33) 

where u is the skew symmetric matrix: 


LU = 


0 

“W 3 

U2 

U> 3 

0 

-wi 

— OJ2 

Wi 

0 


78 


l 



NOTE: 

The angular velocity has to be expressed in the moving reference frame coordinates. 


The kinematical matrix differential Equation 4.33 actually consists of nine first 
order linear differential equations for the nine direction cosines. 

Visualization of Dir ection Cosines 

To provide a visual aid for the rotational motion of a rigid body relative to an 
inertial reference frame, it is convenient to plot the projection of any suitable unit 
vector of the moving body frame on the three orthogonal planes which make up the 
inertial reference frame. 


Z 



Z-X plane: Plot T i3 vs. An 
Z-Y plane: Plot T i3 vs. T ,2 
Y-X plane: Plot .4,2 vs. .4^ 
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a) Orientation Angles 


It was seen in the preceding paragraph that only three independent parameters 
are required for specifying the orientation of a reference frame relative to another 
frame. The scheme to use three independent parameters to define the orientation 
of a reference frame is due to Leonard Euler (1707-1783). It consists of a specified 
sequence of three rotations about three noncoplanar (and nonorthogonal) axes. 
There are two basic types of rotation sequences. 

I) Classical Euler Angles (Type 3-1-3) 



Y 


1) Precession tfr about 3-axis 


' e 'l ‘ 


cos^ sin yj 0 


er 


ei ' 

e 2 

. e 3 . 

— 

— sin ip cos jp 0 
0 0 1 _ 


e 2 

. e 3 . 

= -4 3 (t£>) 

e 2 

. e 3 . 


NOTE: 


There are actually 12 possible Euler rotation sequences considering only positive 
rotations. 

2) Nutation 8 about l'-axis 


e" 1 


‘ l 

0 

0 ' 


' e' ' 


' e 'i ' 

e" 

— 

0 

cos 8 

sin 6 


e' 

= 4 ,( 8 ) 

e 2 

e" 

J 


0 

— sin 8 

cos 0 


. e 3 . 


. . 
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cos (f> sin <f> 0 e" e" 

- sin <j> cos <f> 0 e? = A 3 (<f>) e '2 

0 0 1 J [ e" J e" . 

The matrices A 3 (ip), Ai(0) and A 3 (<f>) are also known as canonical rotation matrices. 
The rotation matrix for the combined sequence of rotations is given by: 

A = A 3 ((f>) Ai{6) A 3 {4>) 

Multiplication of the three canonical matrices yields then: 

(cp exp - scf>c8 sip) (c<j>sip + S(f>c8cip) ( s<ps8 ) 

A = (— s<f> exp — c<j>c8 sip) (—sipsxp + ccpc8cip ) ( c(j>s8 ) 

( s8 sip) (— s6 cip) ( c8 ) 

where c<p = cos^» s<f> = sin<^ etc. 

NOTE; 

The notation angle 8 is usually restricted to the range 0 < 8 < n. 

As in the preceding paragraph it is again possible to establish a kinematical differ- 
ential equation for the rate of change of the Euler angles. To this end we add the 
Euler angle rates vectorically. 

w = i> + e + i> ( 4 . 34 ) 

or in terms of the appropriate unit vectors: 

uii e" 1 + u> 2 e e '2 + u> 3 e'i' = xp e 3 + 8e\ + (f> e$ ( 4 . 35 ) 

The unit vectors on the right hand side can be easily converted by using the above 
canonical rotation matrices. We obtain then the desired relationship in matrix form 
as: 


' 0 ■ 

‘ 8 ‘ 

r 0 ‘ 


= A 3 {4>) Ai{8) A 3 {ip) 0 + A 3 {4>) A\ (0) 

0 

+ a 3 (<p) 0 

(4.36) 

. i> . 

. 0 . 

L 4> . 
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or in component form: 


Wi = ip sin 9 sin <f> + 9 cos <p 
uj 2 = t/> sin 9 cos <j> — 9 sin <p 
lu 3 = ip cos 9 + <f> 

These equations are known as Eulers kinematical equations. 


To obtain the desired differentia] equations for the Euler angle rates the above 
equations are inverted and arranged in matrix form which yields: 


V 

1 

sin <j) 

COS <f) 

0 ' 


’ ^1 

9 

sin 8 

(sin 8 cos <f>) 

(— sin 8 sin (f>) 

0 


U^2 

4 > . 

(— cos 8 sin 4>) 

(— cos 8 cos <f>) 

sin 9 


. W 3 . 


(4.37) 


It is seen that for 9 — 0 the above equations become singular. This singularity has 
its physical origin from the fact that for 9 = 0 the angular velocities <f> and ip (spin 
and precession rates) can no longer be distinguished. This situation corresponds to 
a condition of double gimbal gyroscopic suspensions known as gimbal lock . It turns 
out that the gimbal lock singularity is an intrinsic property of all three-parameter 
orientation schemes. 

NOTE: 


For double-gimbal gyroscopic suspension systems (Cardanic suspension) the nuta- 
tion angle 9 corresponds to the inner gimbal angle whereas the precession angle ip 
corresponds to the outer gimbal angle. 

P.S. A variety of Euler angle systems are in common use. Although there is no 
essentia] difference between them their end formulas (rotation matrix and kinemat- 
ical equation) are difficult to compare. Therefore, it is better to derive each system 
from scratch if needed. 
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II. Modern Euler Angles (type 3-2-1) 



1) Yaw ijj about 3-axis 


‘ e i " 


COS Xp 

sin xp 

0 ' 


ei 


ei 

e' 

= 

— sin ip 

cos xp 

0 


e 2 

= AM 

e 2 

. e 3 . 


0 

0 

1 . 


. e 3 


. e 3 . 


1) Pitch 9 about 2'-axis 



e? ] 


cos 

9 0 - 

- sin 9 


ei ' 


e', ' 


e " 

e 2 

= 

0 

1 

0 


e 2 

= A 7 (9) 

ei 


L C 3 . 


sin 

6 0 

cos 9 


ei . 


e 3 . 

1) Roll <p about l"-axis 








e'/' 1 


' 1 

0 

0 


< ' 


^ 1 


e'" 

e 2 

= 

0 

cos <f> 

sin <f> 


e 2 

= AM) 

e 2 


p'" 

. e 3 J 


0 

— sin <f> 

cos (j) 


e" 

L e 3 . 


p" 

L e 3 J 


The rotation matrix for the combined sequence of rotations is given bj 
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(4.38) 


A = M*) M») Mv) 

Multiplication of the three canonical matrices yields then: 


NOTE: 


.4 = 


(cO cip) ( c9sip ) (sO) ' 

(s<p s9cip - c<j) sip) (s<ps9sip + ctpcxp) (stpcff) 
(c<p s9 cip + sipsip ) (c<p sO sip — s<pcip ) [c<pc9) 


The pitch angle 6 is usually restricted to the range — | < 9 < 

The yaw angle ip is also known as leading or azimuth angle. 

The pitch angle 9 is also known as attitude or elevation angle. 

The roll angle <p is also known as bank or clock angle. 

To obtain the kinematic differential equations for the Euler angle rates we again 
add the angular velocities vectorially: 


u> = ip + 6 + 4>= ip e 3 + 9 e' 7 + <f>e" (4.39) 

By taking the appropriate scalar products as above in Equation 4.35 we obtain the 
desired kinematical equations: 


u>\ = (p — ip sin 9 


^2 = 9 cos <p + ip cos 9 sin <p 


u> 3 = ip cos 9 cos <p — 9 sin <p 


The Euler angle rates are again obtained by inverting the above equation and ar- 
ranging in matrix form with the result: 


' %p ' 

1 

9 


. . 

cos 9 

. 


0 

0 

cos# 


sin <p 

(cos <p cos 9) 
(sin <^>sin0) 


cos <p 

(— sin 0cos0) 
(cos <p sin 9) 



■ U)\ ' 


u; 2 


. ^3 . 


(4.40) 


It is seen that the intrinsic singularity now occurs for 9 = 90°. 
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NOTE: 

The names given to these orientation angles indicate their preferred applications. 
The classical Euler angles are generally chosen for gyroscopic (spinning) systems 
whereas the modern Euler angles find application in the flight dynamics of airplanes 
and missiles. There the interest lies often in small angle deviations (perturbations) 
about a nominal flight condition. It is easily seen that the classical Euler angles 
cannot be used for such situations because of the singularity at 9 = 0. In airplanes 
the body axes are selected such, that the x-axis points “forward”, the y-axis is ”to 
the right” and the z-axis “downwards.” 

C) Quaternions (1843) 

It is possible to avoid the singularity of the three-parameter schemes by spec- 
ifying the orientation of a reference frame using four parameters. This was first 
discovered bv Euler in 1776. These four parameters are therefore also known as 
Euler parameters. When William R. Hamilton (1805-1865) formulated his quater- 
nion algebra it turned out that when applied to rotational kinematics, the quater- 
nions are essentially equivalent to the Euler parameters. In fact, the quaternion 
method proves to be a very convenient tool in the study of rotational dynamics. 

The four-parameter scheme centers about Euler’s theorem: Every rotation of a 
rigid body is equivalent to a single rotation about some axis e (eigen axis) through 
some angle 6. 

It is obvious that any vector lying along this axis of rotation is unaffected by the 
rotation, i.e., it must have the same components in the new and old system. Lsing 
matrix notation for this vector Euler’s theorem can be formulated as: 

R' = .4R=R or(T — £)R = 0 (E41) 

where E = unit matrix and -4 = rotation matrix. 

Recognizing Equation 4.41 as an eigenvalue problem, Euler’s theorem can be re- 
stated as follows: 

The rotational matrix A defining the orientation of a reference frame has the eigen- 
value + 1. 


85 


To prove Euler’s theorem it is therefore necessary and sufficient to show that the 
coefficient determinant 


\A-E |=0 (4.42) 

This can be done by using the orthogonality condition of matrix A and the fact that 
the rotational matrix has the determinant +1. The latter fact can be deduced by 
observing that any rotation matrix must evolve continuously from the unit matrix 
assuming that the two reference frame are coincident prior to rotation. The unit 
matrix E has, of course, the determinant +1. A sudden change in sign from +1 to 
-1 would be incompatible with a continuous motion and represent a transition from 
a right-handed to a left-handed system. 

To arrive at the four-parameter representation we take a geometric approach. Con- 
sider the rotation of a vector T about an axis e through the angle 6 according to 
Euler’s theorem. 



The vector R 2 can be expressed 
in terms of R t and the unit vector 
e along the eigenaxis as: 

R 2 = (Ri • e)e 
+ [Ri — (Ri • e)e] cos S 
+(e x R t ) sin 6 
R2 = Cj + C 2 + C3 


The first term is obviously the projection of along the axis of rotation e. 

The second and third terms have the same magnitude, i.e., Rj sin (R,,e) which is 
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in fact, the radius of the circle which is traversed by the tip of the vector Ri during 
its rotation. 

The second vector lies in the plane of Ri and e, and is perpendicular to e. The 
third vector is perpendicular to the plane spanned by Ri and e. 

It is now possible to convert the above equations to vector-dyadic form by intro- 
ducing the rotation dyadic: 

V = E cos 6 + (ee)(l — cos 6 ) + (e x E) sin 6 (4.43) 

which results in 

R 2 = I> R L (4 44) 


Notice that the rotation dyadic V is defined by four parameters, namely the three 
direction cosines of the rotation axis unit vector e = cosaei + cos/?e 2 + cos 763 and 
the rotation angle 6 . 

To obtain the desired rotation matrix A from the rotation dyadic T>, we have to 
remember that the rotation matrix refers to the orientation of one reference frame 
relative to another. That means, the vector R is assumed to remain fixed and the 
changes in the components of R are due to this rotation of the reference frame. In 
the preceding discussion, it was assumed that the reference frame was fixed and that 
the vector R was rotated in the opposite direction. The rotation matrix, is therefore, 
obtained by simply changing the sign of the rotation angle 6 in Equation 4.43. It is 
important to notice that the direction cosines of the rotation axis refer, of course, 
to the rotating reference frame when they appear in the rotation matrix A. Their 
instantaneous relation to the angular velocity is: 


The result is: 


cos a = 


w 1 
us 


q "2 

cos p = — 
U) 


cos 7 = 


0^3 

U 


A — E cos 8 + ee r (l — cos 6 ) — e sin 6 (4.45) 

where 



cos a 


0 

— cos 7 

cos/3 

e = 

cos /3 

and e = 

cos 7 

0 

— cos a 


cos 7 


—cos ft 

cos a 

0 
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From Equation 4.45 we obtain the elements of the rotation matrix as follows: 


An = cos 6 + cos 2 a( 1 — cos 6) 

(4.46) 

^4 12 = cosacos/?(l 

— cos 6 ) + cos 7 sin 8 

(4.47) 

4 i 3 = cos a cos ^(1 

— cos 6 ) — cos 0 sin 8 

(4.48) 


42 i = COS a cos/?(l — cos £) — cos 7 sin 8 

(4.49) 

4 2 2 = cos 8 + cos 2 0(1 — cos 6 ) 

(4.50) 

4 2 3 = cos/? cos 7(1 — cos 6 ) + cos a sin 8 

(4.51) 

4 31 = cos a cos 7(1 - cos 6 ) 4- cos 0 sin 8 

(4.52) 

432 = cos 0 cos 7(1 — cos 8) — cos a sin 8 

(4.53) 

433 = cos 8 + cos 2 7(1 — cos 8) 

(4.54) 

These elements can be written in simpler form by introducing the Euler parameters: 

. 8 

qi = cos a sin - 

(4.55) 

. . 8 

q 2 — cos p sin - 

(4.56) 

. 8 

q 3 = cos 7 sin - 

(4.57) 

8 

q i = cos - 

(4.58) 

0 i + 02 + 03 + 04 = 1 

(4.59) 

— 7T < 6 < 7T 

(4.60) 

Recalling the trigonometric half-angle formulas 

1 — cos 8 = 2 sin 2 — 
2 

(4.61) 

. ... 8 8 

sin 8 = 2 sin - cos - 
2 2 

(4.62) 
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the first element An can be written in terms of Euler parameters as. 


cos 6 4- cos 2 o(l - cos 6) = 1 - 2 sin 2 ^ + 2 cos 2 a sin ^ 
= 1 - 2(1 - ql) + 2q{ = 1 - 2 (g 2 + q\ + q\) + 2q\ 

= 1 - 2ql - 2 ql = q\ - ? 2 2 - gf + 7? 

The other elements can be written in similar form. 

The notation matrix in terms of Euler parameters is then: 


A = 


(7? - 72 - + 9?) 

2(71 72 - 7374) 
2(71 73 + 7274) 


2(7172 + 7374) 
(-7? + 72 ~ 9a + Qi) 
2(7273 — 7 i 74 ) 


2(7173-7274) 
2(7273 + 7174) 

(-7? - ?2 + ?3 + 94) . 


(4.63) 


NOTE: 

The angle of rotation 6 can be found easily by taking the trace of the rotation 
matrix A in the form given in Equations 4.46 thru 4.54: 


tr A = -4ji + . 4 2 2 4“ ^33 — 1 + 2 cos 8 (4.64) 

It should be noted, that as expected, the four Euler parameters are not independent, 
because one needs only three parameters to define the orientation. They are related 
by the one constraint equation q\ + q\ + 93 + q\ = 1- 

It remains to derive the kinematical differential equations for the rate of change of 
the Euler parameters. Going back to Equation 4.33 we had 


A = —uiA (4.65) 

This matrix equation consists of nine differential equations of which only three are 
independent. For the kinematical differential equations we need four relationships. 
Taking the time derivatives of the diagonal elements of Equation 4.65 and the 
constraint equation associated with the Euler parameters we obtain: 

7i9i - 9292 ~ 7373 + 7494 = ^3(71 92 ~ 9374) ~ ^2(71 73 + 7274) 

- 7 l 7 i + 7272 ~ 7373 + 7474 = -<*>3(71 92 + 9394) + ‘■*■'1(7273 - 7 i 74) 

- 7i7i - 7272 + 7373 + 9494 = ‘*'2(7193 - 9274) - ^1(7293 + 7 i 74) 
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(4.66) 

(4.67) 

(4.68) 


01 01 + 0202 + 0303 + 0404 = 0 


(4.69) 


Introducing the matrices 


■ 01 ' 


“ 0 


U>2 

“^1 



02 

UJ — 


0 

-U>i 

— UJ2 


Cj 1 — UJ 

03 


— 

CUi 

0 

~^3 



. 04 . 


. ^1 

UJ2 

4J3 

0 


. w | 0 


the above equations can be written in matrix form as: 

1 _ 

q = (4.70) 

Integration of Quaternions 

The kinematical Equations 4.70 are not numerically integrated as such. A better 
integration procedure can be developed by assuming for a moment that the angular 
velocity us is constant. For this case the quaternions of Equations 4.55 thru 4.58 
can be integrated in closed form, because the direction cosines remain constant, 
such that at all times: 

(4.71) 

angle 6 = u> t we 


U>2 

cos a = cos 3 = — cos 7 = 

UJ u 


w 3 


Setting the initial rotation angle 6q — u t$ and the final rotation 
have then: 


^1 . 6 + 6 q 

?i = — sin — - — 
uj 2 

^2 . 6 + 

92 = — Sin — — 
uj 2 

^3 . 6 + 6 0 

9 a = — sin — - — 
uj 2 


9 4 = cos 


<5 + by 


0 


90 



To arrive at the relationship between the initial and final quaternions, let us also 
introduce: 


wi . 

9io = — sm 

id 

id? 

920 — sin 

id 

930 — sin 

id 


2 

So 

2 

So 

2 


940 = cos 


S _ o 
2 


The next step will be illustrated using q x as an example. Expanding the trigono- 
metric function, we obtain: 


9 1 = 


uh 

id 



6 ° 

cos — + cos 



id\ . 8 J So, S .id 1 . 6 o\ 
91 = — sm - (cos -) + cos - (-sin —) 


id 


id 


id i . S So 

q i - — sin - 940 + cos —910 
id l £ 

The other components are manipulated in like manner. We obtain them: 

S ui . S 

91 = cos 2 910 + Sin 2 940 

5 W 2 . S 

92 = cos - 920 H sin 7, 940 

L id l 

6 ids . S 

q 3 = COS - 9 30 + — s,n o 940 

/ LJ / 

6 , 0 J\ 0J2 ^ 

o 4 = — sin — ( — 910 H 920 + 930) + cos 940 

’ 2 w u> w l 

where 6 = ut, 


or in matrix form: 


. u>t u.ui 

q = (E cos — - - sin — )q 0 
2 u l 
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where £ = (4x4) unit matrix and 57 defined as in Equation 4.70. While only an 
approximation for a time-varying angular velocity the error is extremely small for 
short time integration intervals. This scheme automatically satisfies the constraint 
equation and its computer efficiency is excellent. 

NOTE: 

Computer coding is more efficient when using the above components because the 
matrix equation involves some cancellation for the off-diagonal terms. 

d) Interrelatedness 

Quite often it is necessary to convert from one orientation scheme to another. 
The conversion from Euler angles or quaternions to direction cosines is, of course, 
directly given by the rotation matrix .4 expressed either in Euler parameters or 
Euler angles. The inverse problem of converting from the direction cosine (rotation) 
matrix to Euler angles or quaternions is more involved due to the inherent ambiguity 
resulting from having more equations than unknowns. 

A) Classical Euler Angles from Rotation Matrix 

9 = COS -1 .4 33 0 < 9 < X 

<t> — tan 1 A13/A23 — 7T < <f> < 7r 

V* = tan 1 A 31 /A 32 — 1 r < < x 

B) Modern Euler Angles from Rotation Matrix 

9 = sin -1 (-.4 33 ) — tt/2 < 9 < x/2 

<f> = tan 1 A 23 /A 33 — x < <f> < x 

t/j = tan 1 A12/A11 —x < \p < x 


The above relationships can be easily verified by inspection of the rotation matrices 
expressed in terms of Euler angles. The ambiguities stemming from the inverse 
tangent functions can be resolved by observing the signs of the direction cosines in 
their arguments as numerators and denominators. 

The following quadrants can then be determined: 

I. Quadrant: +/+ 
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II. Quadrant: +/- 

III. Quadrant: -/- 

IV. Quadrant: -/+ 




6 ip 6 . ip . 6 <f> 

q i = sin - cos - cos - - sin - sin - cos - 

0 ip <f> . fp . 8 

q 7 = sin - cos — cos - + sin — sin - cos - 

ip 0 ip . 0 . <f> tp 

q 3 = sin — cos - cos - - sin - sin - cos — 

tP 6 . <p ip 6 <p 
q 4 = sin - sin - sin - + cos - cos - cos - 
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KQTE; 

For small angles the quaternions are approximately: 


_ <t> 
qi 2 

9 

q 2 = - 

2 

_ ij) 

93 — 2 

94 = 1 

Notice also the ubiquitous presence of half angles which is typical for all four- 
parameter schemes. 

The derivation of these formulas can be done taking the following steps: (Quaternion 
algebra would provide a more elegant way). 

Take the diagonal terms of the rotation matrix .4 and the constraint equation of 
the quaternions: 

9 ? - ql - 93 + 94 2 = - 4 „ 

~ 9l "f ?2 ~ ?3 F ?4 = -422 

-q] - 92 + 93 + 94 = -433 

ql + ql + 93 + 9 4 2 = 1 


As an example we solve for q 7 by adding the second and fourth equation to obtain: 

2 ^2 2 9 4 — 1 + A 22 

Also adding the first and third equation: 

(4.72) 

— 2 + 2 ql ~ An + A 33 

Eliminating q\ yields then: 

(4.73) 

Aq 7 = 1 — - 4 n + -4 22 — -4 33 

Next we introduce the half-angle relations: 

sin q 2 = 2 sin — cos — 
2 2 

(4.74) 
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COS £*1=1—2 


sin 


<*i 

2 


To make the notation more concise we let sin c*i = «i, costt] = Cx,sin = s t and 


cos ^ = C\ . 


Introducing in addition, the notation for the modern Euler angles: xp = 0 U 8 = 9 2 
and <f> = 6 3 , we rewrite the direction cosines as: 


An = ci c 2 

A 22 = Si S 2 S3 + Ci C3 

.433 = C 2 C 3 


The second quaternion can then be manipulated as follows: 


4 q\ 

— 

1 + 

s 2 ^3 

+ Cl c 3 - 

•Cl c 2 - c 2 c 3 



— 

1 + 

8 S 1 S 2 S3 Ci c 2 c 3 

+ (1 -2«?)(1 - 2 s 

3) 


- 

(i- 


-2*1)- 

(1 — 2 s 2 )(l — ^ S3) 



= 

2 Si 

s 2 S3 Cl 

c 2 C3 + S 2 

+ sf S3 — s 2 si — s\ 

4 

= 

2 si 

S 2 S3 Cl 

c 2 c 3 -1- Sj 

{s\ + c 2 ) (s 2 + c 2 ) 



+ 

$1 s 

*(•§+< 

*) - , 2 s 2 ( 

+ 4 ) - + 

4 ) 


— 

2si 

S 7 S3 Cl 

C 2 C3 + S2 

c 2 cl -1- s\ S3 C^ 



Finally: 


Taking the square root we 


q\ = (s 2 C] c 3 + Si s 3 c 2 ) 2 (4-75) 

take the positive sign and that for small angles q 2 = 


<72 — $2 C 3 + s l s 3 c 2 
or reverting back to the original notation: 

6 xj) (f> . xj) . <f> 0 

q 7 = sin - cos - cos - + sin — sin - cos - 

The other quaternions can be obtained in a similar fashion. 


(4.76) 


(4.77) 


E) Quaternions from Rotation Matrix 

Different algorithms have been suggested to solve this inverse problem. The most 
elegant an efficient one makes use of the inherent symmetry 

of the quaternion relationships. To reveal this symmetry we introduce the fol- 
lowing definitions: 


p = 2 q, An = 7/t and Tr = A n + A 22 + ^33 (4-78) 
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Weakness: 


The kinematical equations are transcendental and singular (gimbal look). As a 
consequence, computer coding is ineffective and awkward. 

Quaternions: 

Strength: 

The small set of four parameters having linear nonsingular kinematical equa- 
tions admits of highly efficient computer coding. Quaternion algebra affords concise 
derivation of multiple reference frame inter-relationships. 

Weakness: 

Poor visualization. 

E) Direction (Aerodynamic) Angl es 

It was obsened that the definition of an orientation (attitude) requires three 
independent parameters. By contrast a direction can be specified using only two 
independent parameters. This can be understood by realizing that a direction is 
defined by a unit vector whose three direction cosines have to satisfy the constraint 
equation cos a + cos 2 3 + cos 2 7 = 1. This leaves two independent parameters 
for the direction specification. Direction angles find their foremost application in 
determining the aerodynamic forces acting on a body. This explains the terminology 
commonly used for the direction angles. There are two types of direction angles in 
common practice. It is assumed that prior to the rotation sequence the x-axis of 
the body is aligned with the desired direction. In aerodynamic terms, the x-axis is 
aligned with the nominal flight velocity v. 

Type I: 


Here the direction angles are defined by the following rotation sequence: 
1) Counterclockwise rotation /? about z-axis 


[e; 1 


cos 3 

— sin 3 

0 

*2 

= 

sin 3 

cos 3 

0 

. ej J 


0 

0 

1 


2) Clock wise rotation a about y '- axis 
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' cos a 0 -sin a 1 


’ e i 


O 

o 


e 2 


sin a 0 cos a 


. e 3 . 


The total rotation matrix is obtained by multiplying the two canonical rotation 
matrices: 

e" " cos a cos 8 — cos a sin 8 — sin a 
e" = sin 0 cos 0 0 

ej sin a cos 0 — sin a sin 3 cos a _ 

As assumed the flight velocity before the rotation was along the x-axis: (v = v\) 
Therefore: 


i 

j 

k 


v = 


V 

0 

0 


The velocity components in the rotated reference frame are therefore. 


v x = v cos a cos 3 
v y = v sin 3 
v z = v sin a cos 3 


The two orientation angles a and 3 can now be expressed in terms of these velocity 
components as: 

Angle-of-attack tana = ^ -180° < « < 180° 

Sideslip angle sin 3 — if < 3 < 90° 


The restrictions on the ranges of the orientation angles are introduced to avoid 
ambiguities when taking the inverse trigonometric functions to obtain the direction 

angles. The quadrant of the angle of attack 

is determined by the appropriate signs in the numerator and denominator of the 

argument. 

Comparing the direction angles with the modern Euler angles, it is readily seen 
that the angle-of-attack corresponds to the pitch angle 9 and the sideslip angle 0 to 
the negative (!) yaw angle V>. A positive sideslip angle 3 obtains for a nose left 
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situation. This unusual sign convention stems from the desire to have a positive 
sideslip angle when the wind comes from the right . 

Type II: 


Here the direction angles are defined by a slightly different rotation sequence which 
essentially corresponds to the modern Euler angle rotation sequence with = 0. 

1) Clockwise rotation a T about y-axis 





' e i ' 


cos ax 0 - 

- sin a T 

i 




e' 

— 

0 1 

0 

j 




. e 3 . 


sin a T 0 

cos a R 

k 

2) Clockwise 

rotation 4> R about 

I'-axis 





r 

e" 1 


' 1 0 

0 ' 

e > 1 




e" 

= 

0 cos <f> R 

sin <p R 

e' 



* 

e 3 J 


0 — sin <p R 

COS (f) R 

. e 3 . 

The total rotation 

matrix is 

then: 






■ 

cos ax 0 

— sin aj 




sin 

ax sin <j) R cos <j> R 

cos q x sin (f> R 


i p" 

L J 


[ sin otT<t>R —s\r\<j) R cosa T cos<j) R 


Assuming again that the nominal flight velocity is along the original x-axis. the 
velocity components in the rotated frame are: 

v x = v cos a T 
v y = v sin aj sin (f> R 
= v sin a T cos <j> R 

The two orientation angles can then be defined as: 

Total angle-of-attack cosqj- = ^ 0 < aj < 180° 

Roll Angle tan <p R = 0 < <f> R < 360° 

Notice again the range restriction for the total angle of attack which allows a 
unique inversion of the cosine function. 
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cos Qj = 1 — 2 sin 2 


2 


To make the notation more concise we let sin ai = s 1; cosai = ci,sin ^ = s y and 

COS ^- = Ci . 


Introducing in addition, the notation for the modern Euler angles: = 0i,0 = 

and <j> = 8 3 , we rewrite the direction cosines as: 

Tn = C 1 c 2 
A 22 = 5] S 2 S 3 + Ci C3 

-^33 = C 2 C 3 

The second quaternion can then be manipulated as follows: 


4 < 7 2 

= 

1 + 

$1 $2 

+ C 1 C3 - 

- Cl c 2 - 

c 2 C 3 


= 

1 + 

8 s { s 2 S3 Ci c 2 c 3 

+ (1-2 

^)(l-2s 2 ) 


— 

(1- 

-2s 2 )(l 

-2 s 2 )- 

(1— 2 s 2 

)(l-2s 2 ) 

92 

= 

2 si 

S 2 S3 Cl 

c 2 c 3 + si 

+ s?s 2 - 

2 2 2 2 
- S\ S 7 — Sj S3 


= 

2 si 

s 2 s 3 C 1 

C 2 C 3 + S2 

W + <$) 

( S 3 + C 3 ) 


+ 

s 2 s 

2 (s 2 + r 

: 2 ) - s 2 s 2 ( 

si + cl) - 

- si«l(s 2 + c 2 ) 


= 

2sj 

s 2 S3 Cl 

c 2 C 3 - 1 - s\ 

c 2 cl + s 2 

s§ c\ 


Finally: 

q\ = (s 2 Ci c 3 + Si s 3 c 2 ) 2 (4.75) 

Taking the square root we take the positive sign and that for small angles q 2 = |: 

92 = S 2 Ci C 3 + Si s 3 c 2 
or reverting back to the original notation: 

8 \j) <j> . ip . <j> 8 

o 2 = sin - cos — cos — + sin — sin — cos - 

H 2 2 2 2 2 2 

The other quaternions can be obtained in a similar fashion. 


(1.76) 

(4.77) 


E) Quaternions from Rotation Matrix 

Different algorithms have been suggested to solve this inverse problem. The most 
elegant an efficient one makes use of the inherent symmetry 

of the quaternion relationships. To reveal this symmetry we introduce the fol- 
lowing definitions: 


p — 2 q, A 4 4 — Tr and Tr — An + T 22 + T 33 (4.78) 
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where Tr is the trace of the rotation matrix. The diagonal terms of the rotation 
matrix Equation 4.63 can be cast into a set of symmetrical equations: 


Pi = 1+2 A n - T r 
p\ = 1 + 2 .4 22 — Tr 
p\ — 1 + 2 .4 33 — Tfi 
p\ = 1 + 2 .4 44 - T r 


(4.79) 


The off-diagonal terms of Equation 4.63 furnish all combinations of cross product 
terms: 


Pl P2 — *4i2 + - 

4 2 i 

Pi Pi = • 

4 12 

“ A 21 


Pi Pi — V 3 i + .■ 

T 1 3 

P 2 Pi = ■ 

431 

“ -4 13 

(4.80) 

P‘2 Pi = *4 23 + - 

^32 

Pi Pi = - 

423 

— A 32 



The first step in solving the above equations for the unknown quaternions is to 
select the largest p, of the diagonal Equations 4.80. 

This assures the highest possible numerical accuracy for the solutions. The next 
step is to choose the appropriate three off-diagonal equations from Equation 4.80 
to calculate the remaining three unknown quaternions. 

To avoid ambiguities the positive square root is taken when the first quaternion is 
calculated in step one. It is customary to restrict the range of the rotation angle 
to —7 r < 8 < +7r. From the definition of the Euler parameters in Equation 4.58 it 
follows that q A is always positive. This requires that all the signs of the quaternions 
have to be changed if the above algorithm yields a negative p 4 . A useful computer 
code for the solution of the above problem can be set up by taking the following 
steps: 

Step 1: Define matrix 


An 

A 12 

^13 

A 23 

A 2 l 

A 22 

*4 23 

A 3 1 

A 31 

A 32 

*4 33 

A 12 

CN 

CO 

f 

— j 

”A i3 

— *421 

1 
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Step 2: Define matrix 


p = p + p T + (1 - T r )E 
where E is a 4 x 4 unit matrix. 

Step 3: Select k = MAX p„ 

where k is the row of the largest diagonal element of p. 
Step 4: q] = Pkj/^^/Pkk 
Step 5: qj = q] sgn q‘ 
where sgn = — 1 for qj < 0 
+ 1 for ql > 0 


Summary: 


Direction Cosines: 


Strength: 

Linear nonsingular kinematical equations. Easy visualization makes them well 
suited for attitude display plots. 

Weakness: 

The unduly redundant set of nine parameters requires excessive computational 
effort. 

Euler Angles: 

Strength: 

The set of three independent parameters makes them a natural tool for analytical 
studies. The classical system is particularly useful for analyzing 

gyroscopic (spinning) systems whereas the modern system is applied to the sta- 
bility and response analyses of systems which deviate only moderately from nominal 
operating conditions. For these cases the classical system would become singular. 
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Weakness: 


The kinematical equations are transcendental and singular (gimbal look). As a 
consequence, computer coding is ineffective and awkward. 

Quaternions: 

Strength: 

The small set of four parameters having linear nonsingular kinematical equa- 
tions admits of highly efficient computer coding. Quaternion algebra affords concise 
derivation of multiple reference frame inter-relationships. 

Weakness: 

Poor visualization. 

E) Direction (Aerodynamic) Angles 

It was observed that the definition of an orientation (attitude) requires three 
independent parameters. By contrast a direction can be specified using only two 
independent parameters. This can be understood by realizing that a direction is 
defined by a unit vector whose three direction cosines have to satisfy the constraint 
equation cos 2 a + cos 2 5 + cos 2 7 = 1. This leaves two independent parameters 
for the direction specification. Direction angles find their foremost application in 
determining the aerodynamic forces acting on a body. This explains the terminology 
commonly used for the direction angles. There are two types of direction angles in 
common practice. It is assumed that prior to the rotation sequence the x-axis of 
the body is aligned with the desired direction. In aerodynamic terms, the x-axis is 
aligned with the nominal flight velocity v. 

Type I: 


Here the direction angles are defined by the following rotation sequence: 
1) Counterclockwise rotation 0 about z-axis 


'-i 

*2 

ei 


COS ft 

— sin d 

0 ' 


i 

sin ft 

cos 3 

0 


j 

0 

0 

1 _ 


k 


2) Clockwise rotation a about y'-axis 
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It is easily observed that the type II orientation angles correspond to the last two 
rotations of the modern Euler angle system. Therefore the total angle of attack a T 
can be identified with the pitch angle 6. Whereas the roll angle <j> R is really nothing 
else but the role angle <f>. 

Both systems of orientation angles appear to be equally in use. However, the type 
II angles are not suitable for dynamic studies of small perturbations from a nominal 
flight condition because the roll angle approaches an indeterminate form (v y = v z « 
0). It can be seen that the type 1 angles are actually well behaved for this condition. 
Notice also that in this case the x-component of the velocity v x is nearly equal to 
the nominal flight velocity v (v x ss v). 

The types of orientation angles are mathematically related as follows: 


A) Conversion from type I to type II: 


cos a R = cos a cos 0 
ig<i>R = (sin a)~ 1 tg0 


B) Conversion from type II to type I: 


tga = cos <f) R tgaT 
sin 0 = sin < j> R sin a? 


Extreme caution has to be exercised when converting from one type to the other to 
make sure that the orientation angles lie in their correct quadrant. 

For large orientation angles the type II system is easier to visualize than the type I 
system . 

NOTE: 


If the orientation angles have to be determined by starting from a given direction, 
the rotation sequence has to be reversed and the rotations be performed in the 
opposite clock sense. Thus the type I system would require first a counterclockwise 
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rotation about the y-axis and then a clockwise rotation about the z-axis to align 
the given direction with the positive x-axis. Using the type II system, the reverse 
rotation sequence starts with a counter clockwise rotation (roll) about the x-axis 
followed by a counter clockwise rotation (pitch) about the y'-axis. 
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4.4 Moment of Inertia Properties: 

The following properties of the moment of inertia matrix and its elements (moments 
and products of inertia) play a central role in dynamics studies. 

A) Triangle inequality of the moments of inertia 

The three moments of inertia about three rectangular axes are such that the 
sum of any two is greater than the third. 

Let: 

I n = ^{ X l + * 3 )™ 

7 22 = £(zi + xl ) 171 
7 33 = £(z? + x\)m. 

Then: 


1 11 + 1 22 — 1 33 

= £(*2 + xl 4- x\ + z 3 - x\ - x\)m = 2-Kx\m > 0 
Therefore: 7 n + 7 22 > 7 33 Q.E.D. 

Similar relationships are obtained by cyclic permutations of the indices. 

B) Theorem of Parallel Axes (Steiner Theorem) 

Let O be the origin of the reference frame about which the moment of inertia dyadic 
is defined and C.M. be the mass center of the rigid body. 
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In the figure: 

d = distance of C.M. from O 
r, = location of mass element m,-. 

Using vector-dyadic notation the inertia dyadic about O can be written as: 

Jo = £ [R-£ - (R, R,)] m, 

Jo = E [(d + r ifE - (d + r.) (d + r,)] m, 

Io = E [(d 2 -I- 2d • r, 4- r 2 )£ — (d d 4- dr, + r,d -I- r, r,] m, 

For the mass center Er, m t = 0 and Ed m l = M d where M = E m,. 

Therefore: 


I 0 = E(r 2 <? - r, r;) m, + M(d 2 £ - d d) 


or 

lo = Ic.M. + M{d 2 S - d d 


Converting to matrix form, we introduce: 


d — d x e x + c? 2 e 2 + ^ 3 e 3 

It can be easily verified that the dyadic C = a 2 J — aa corresponds to the matrix 

L = —a 2 

where 


a = 


0 

-as 

a 2 

<*3 

0 

-a\ 

“02 

a i 

0 


Therefore the matrix form of the above inertia dyadic T 0 is: 


h ~ I c.m. — M d 2 
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In component form: 


hi = I\ x + M{d\ + d\) 

I22 = 1 22 T hi (d 1 T ^3) 
hi = hi + M(d\ + d\) 

and also: 

I12 — 1 12 ~ M d\d.2 
hi = I‘a-Md id 3 
hi = hi ~~ hi ^2 c/3 

It can be seen that a translation of the axes away from the mass center always 
results in an increase of the moments of inertia. On the other hand, the products 
of inertia may increase or decrease depending on the particular situation. 

C) Theorem of Rotated Axes 

Here we establish the relation between the moment of inertia matrices expressed 
in two different reference frames. To do this we calculated the rotational kinetic 
energy for the two different systems whose common origin is 

chosen at the mass center. For convenience sake we call one system the primed 
system and the other the unprimed system. It is apparent that the rotational kinetic 
energy, being a scalar, has to be the same for both coordinate systems. 


T = ^ u )' T I'u)' = ) T Iu) 

The rotational transformation is: 


u>' = Auj 


The kinetic energy is therefore: 

T = ]-u T (A T I 1 A)u> = ^u) 7 (/)u> 
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The two bracketed terms must be identical, i.e., 


A t I'A = I 

By proper pre-and post-multiplication with the rotation matrix A, we can solve for 
the primed moment of inertia matrix: 


I' = AIA t 

This is the desired transformation law of the inertia matrix when going from an old 
(unprimed) to a new (primed) system by a rotation. 

D) Principal Axes 

The moment of inertia matrix is a real symmetric matrix. It is a theorem of 
matrix analysis that any real symmetric matrix can be reduced to diagonal form 
by means of an orthogonal transformation. This means that we can always find 
a reference frame in which all products of inertia are zero. This reference frame 
is known as principal axes system and the three mutually orthogonal coordinate 
axes as principal axes . The diagonal terms of this diagonal inertia matrix are called 
principal moments of inertia . 

The problem of finding such a principal axes frame is equivalent to solving the 
eigenvalue problem 


(I -\E)x = 0 

The eigenvalues of this problem are the principal moments of inertia 

and their associated eigenvectors represent geometrically (unit) vectors along 
the principal axes of the reference frame. The eigenvectors can be grouped 
together to form a matrix $ whose columns are the eigenvectors x 1 ,x a , and x 3 . 

Designating the components of the eigenvectors as 


^ll 


’ *12 ' 


*13 ' 

X 2I 

X 2 = 

*22 

X 3 = 

*23 

. *31 . 


. *32 . 


. *33 . 


the matrix of the column vectors called modal matrix is 
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= [xi|x 2 |x 3 ] 



x n 

X 12 

*13 

# = 


X 22 

*23 


- ®31 

^32 

*33 . 


Let also A t = 7i, A 2 = 7 2 . ^3 = h- 


The diagonalization of the inertia matrix I using the modal matrix # can then be 
performed as follows: 


7 # 


[7*i 

[7i*i 


*1 


Ix 2 | 7* 3 ] 

72*2 | 7 3 aj 3 ] 


*2 


*3' 


7j 0 0 

0 7 2 0 

0 0 7 3 


= *I D 


The second step in the above equation is a consequence of the eigenvalue solution 
Ixi = Ai Xi = I\ X\. Premultiplying by # T yields: 


7 # = # T # 1 D = Id 


or finally: 


I D = 7 # 

Comparing this transformation law with the law governing the change of the inertia 
matrix under rotation we see that the modal matrix is related to the rotation 
matrix A simply as: 


.4 = <f 7 


The rotation matrix is the transposed modal matrix. This rotation is also known 
as principal axes transformation . 

Example: 

For the Gemini spacecraft the inertia matrix referred to the control axis system 
with origin at the mass center was: 
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I = 


4560.7 31.2 -43.4 ' 

31.2 4545.0 -270.7 

-43.4 -270.7 1567.4 


We want to make a principal axis transformation: 


1) Solve the eigenvector problem 


Eigenvalues: 


(I - \E)x = 0 


/, = 4530.11944 
h = 4600.531811 
h = 1542.448749 


Modal Matrix: 


= 


0.74693 

-0.66305 

0.04964 


0.66477 

0.74319 

-0.07584 


0.01339 

0.08965 

0.99588 


NOTE: 


The eigenvectors are normalized to unit length. 

2) Check the sign of the determinant of #. 


| # |= 0.9999997403 

If | # | turns out to be negative we change the sign on one row or column. 

3) Find minimum rotation angle 6: 


tr $ = 1 + 2 cos 6 

To do this we have to change the sign on any two rows. 


108 



tr #! = +0.74693 + 0.74319 + 0.99588 — 

tr # 2 = -0.74693 - 0.74319 + 0.99588 — 

tr # 3 = -0.74693 + 0.74319 - 0.99588 -> 

tr # 4 = +0.74693 - 0.74319 - 0.99588 -* 


6 = 42.01° 
6 = 138.34° 
6 = 178.88° 
6 = 174.19° 


4) The rotation matrix A is obtained by taking the transpose of the modal matrix 
having the minimum rotation angle 6 : 


A = #f = 


' 0.74693 
0.66477 
0.01339 


-0.66305 

0.74319 

0.08965 


0.04964 

-0.07584 

0.99588 


The orientation of the principal axis system relative to the control axis system can 
be expressed in terms of the classical Euler angles or in terms of the modern Euler 
angles. 


Example: 

Determine the direction of the minimum principal moment of inertia axis e 3 relative 
to the corresponding control axis. We obtain: 


9 = cos 1 j 4 33 = cos 1 0.99588 


9 = 5.2° 


<fi = tan 1 


-4 13 _i 

= tan 

+23 


0.04964 

-0.07584 


This angle lies in the second quadrant. 


<f> = (180 - 33.206)° ^<j>= 146.8° 


E) Ellipsoid of Inertia 

The inertial properties of a rigid body can be conveniently depicted by an ellipsoidal 
surface which is in essence a plot of the moment of inertia of the body as a function 
of the rotation axis direction. Going back to the matrix form for the rotational 
kinetic energy of a rigid body we have 
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T = ^u> T Iu> = ^I A u 7 (4.81) 

where the scalar I A is the moment of inertia about the instantaneous rotation axis. 
The single scalar expression for the kinetic energy on the right side can be verified 
by letting the angular velocity u; coincide momentarily with one of the coordinate 
axes. Now let us define a vector p having the same direction as u such that: 

1 

P = 7f=U 

UJy/I A 

Inserting this in Equation 4.81 above, we obtain 


p T I p — 1 where p = 

Consider now p to be a positive vector drawn from the origin O to a point (x, y, z) 
then we can write the scalar equation: 



hix 2 + / 222/ 2 + Cs3-z 2 + 2/ 12 ar y + 2/ 13 i z + 2/ 2 3j/2 = 1 

This ellipsoidal surface centered about O is called the ellipsoid of inertia . 

The coordinate transformation which brings the inertia ellipsoid in its standard 
form is, of course, exactly the principal axis transformation previously discussed. 

The moment of inertia about any rotation axis can be found directly from the 
magnitude of the vector p. 


I A = 


1 

Tp 


where p is the length of the straight line drawn from the origin 0 to a point on the 
surface of the inertia ellipsoid. 


NOTE: 


A quantity closely related to the moment of inertia is the radius of gyration k A 
defined by 


no 



I A = Mk\ 


In terms of the radius of gyration p can be written as: 

1 

p = 7= UJ 

uk A VM 

The inertia ellipsoid is fixed with the body and rotates with it. Two rigid bodies 
having equal mass can have quite different shapes and still have identical ellipsoids 
of inertia. They are called dynamically equivalent or equimomental. Although the 
inertia ellipsoid is a figure of revolution, the corresponding rigid body need not be. 
Roughly speaking the shape of the inertia ellipsoid is similar to the shape of the 
corresponding rigid body. For instance, a prolate (oblate) rigid body has a prolate 
(oblate) inertia ellipsoid. The standard form of the inertia ellipsoid is 

I\x\ + I 7 x\ + /3Z3 = 1 

where / t , I 2 and / 3 are the principal moments of inertia. 
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4.5 Free Motion of a Rigid Body 


The free motion is characterized by the absence of external moments. The origin 
of the reference frame is either fixed in space or at the mass center. Both cases are 
dynamically identical. There are two ways of treating this problem - the geometric 
(Poinsot method 1834) and the analytic (Euler case 1758) way. 

To simplify the problem, the reference frame is chosen to be a principal axes 
system. There is, of course, no loss of generality involved in this particular choice 
of the reference frame. 

The equations of motion are obtained from Equation 4.12 and read in scalar 
form: 


I\W 1 + (^3 — ^2) w 2^3 — 0 


I ^ UJ *2 + (^1 — ^3)^1 ^3 — 0 


(-1.82) 


^ 3^3 + (^2 — — 0 


A) Poinsot Method (1834) 


This method represents a geometric solution of the torque-free motion of an 
unsymmetrical rigid body. W ith no external torque acting on the body, the kinetic 
energy and the angular momentum are conserved. 

Therefore: 


2 T — + ^ 3-^3 (4.83) 

This equation represents geometrically an ellipsoid known as Poinsot Ellipsoid or 
kinetic energy ellipsoid. It differs from the standard inertia ellipsoid only by the 
scale factor \f¥T . 

Also: 


H 2 


Jfw? + + I’fai = 2 TD 


j 2 2 


2 2 _ 


This ellipsoid is the angular momentum ellipsoid . 


(4.84) 
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The quantity D is defined as: 


D — — = CONSTANT (4.85) 

It has the dimension of a moment of inertia and is used in the subsequent discussion 

as a matter of convenience. 

The Poinsot method does not yield the angular velocity a? as a function of time but 
only the path traced by the instantaneous rotation (spin) axis on the rigid body. 
This trace is called the polhode and is the curve obtained by the intersection of the 
Poinsot ellipsoid and the angular momentum ellipsoid. Mathematically speaking, 
it is given by the simultaneous solution of Equation 4.83 and Equation 4.84. 

For the ensuing discussion, we assume without loss of generality that 

h<h< h (4.86) 

It can then be shown that for this case, the constant D must lie in the range 

h < D < h (4.87) 

As a consequence, a polhode angular momentum ellipsoid is always more elongated 
than the corresponding polhode Poinsot ellipsoid. The following two figures illus- 
trate the geometrical aspects of the Poinsot method. 

It is observed that the polhodes form closed curves about the smallest and largest 
moment of inertia reflecting a stable motion. On the other hand, the polhodes in 
the vicinity of the intermediate moment of inertia / 2 have hyperbolic character 
indicating an unstable motion. 

In the case of axial symmetry, the Poinsot ellipsoid and the angular momentum 
ellipsoid become ellipsoids of revolution (spheroids). The polhodes become circles 
perpendicular to the spin axis. The rotation about the axis of symmetry becomes 
more stable, whereas the rotational motion about a transverse principal axes is 
neutrally stable. 

The polhodes furnish an argument for the stability behavior of the rotational 
motion about a principal axis in the presence of internal energy dissipation. 

Since 


h<{D=—)< h (4.88) 

from Equation 4.87 it is observed that in this case the kinetic energy T must decrease 
while the angular momentum H remains constant. As a consequence, the quantity 
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Figure 4.1: Intersection of Poinsot and Momentum Ellipsoid 
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D is increasing which means that a motion in the vicinity of the principal axis of 
minimum moment of inertia (D ~ 7\), will gradually go over to a rotation about 
the principal axis of maximum moment of inertia (D ~ I 3 ). 

Thus, the principal axis of minimum moment of inertia is one of unstable equi- 
librium in the presence of internal energy dissipation. This behavior was actually 
observed for some satellites notably the Explorer I satellite which was spin-stabilized 
about the longitudinal axis of minimum moment of inertia. The source for the en- 
ergy dissipation can be provided by internal fluid motion or feasible antennas. The 
important question is how long a spinning satellite without large changes in orien- 
tation. 

The Poinsot method can also be used to calculate the limits of the mutational 
or wobble motion of an unsymmetrical torque-free rigid body. To this end, we 
write the Poinsot ellipsoid (Equation 4.83) and the angular momentum ellipsoid 
(Equation 4.84) in terms of the angular momentum components. Thus, 

IJ2 IJ 2 r/2 

-y- + — 2T (4.89) 

1 1 h h 

//j 2 + H% + H% = H 7 Momentum Sphere (4.90) 

Multiplying Equation 4.89 by D and observing that by definition H 7 = 2 TD we 
can combine both equations to obtain the path of the angular momentum vector 
on the momentum sphere: 

JT?(1 - r ) + H\( 1 -y)+ H\( 1 - j ) = 0 (4.91) 

h -*2 -*3 

The trace of the momentum vector H on the momentum sphere is illustrated in the 
figure. 

Case A: Spin about I i-axis ( D < 1 2 ) 

For 9min : H 2 = 0(w 2 = 0) 


Hi _ 1-D/h _ h (D-h) 
Hf 1 - D/h h (h ~ D) 

is h,,N = 

For 6 max '• H3 = 0(013 = 0 ) 
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Figure 4.3: Trace of H on Momentum Sph 


^2 _ 1 ~ D / h _ h , D — 1 1 
H? ~ 1 -D/h ~ h { h~D ] 

is h,ax = 

Spin about / 3 -axis ( D > / 2 ) 

< '■ H% = 0(u2 2 — 0) 


HI. 

HI 


v ■ H 1 — 0(u.'i — 0) 


1 - D/h 
1 - D/I , 

*5 @ MIN = 


— ( h ~ ^ 

h ^D-h 

ll ( H — D 

l*D-h 



) 


) 


Hi 


1 ~ J >// 3 

1 - D/h 

tg $ MAX = 


— ( H 
h 1 t >-/ 2 

li ( h ~ ^ 
/ 3 l z >-/ 2 


) 


) 



B) Analytical Method 


The analytical solution of the Equation 4.82 for the torque-free unsymmetrical 
body was first obtained by Euler in 1758. He showed that the angular velocity 
components become elliptic functions of time. 

There are, of course, also three very simple particular solutions of Equation 4.82 
namely: 

w, = CONSTANT w 2 = = 0 

w 2 = CONSTANT uj\ — — 0 

w 3 = CONSTANT Wi = w 2 = 0 

These solutions represent steady rotations about the principal axes of inertia. These 
are the only axes about which the body will spin steadily. 

The derivation of the Euler case will not be given only the solution. 

The elliptic functions appearing in the solution are defined as follows: 

The elliptical integral of the first kind is: 

r y dy 

U ~Jo ^(l-y^l-fcV) 

that means u is a function of y and k: 


u = F{y, k) 

The elliptic function is then the inverse function 


y = F~ l (u, k) = Sn ( u,k ) 
where k is called the modulus (0 < k < 1) 

The solution for the angular velocity components can then be written as. 


w, = H 


N 


(D - h) 

1 , 0 ( 1 ,-!,) 


Vn (A <, k) 
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w 2 = -H 
^3 = H 




(Ji ~ D) 

hD(l, - / 2 ) 


N 


i!> - £1 
1 , 0 ( 1 , - /,) 


Sn (A/, k ) 


Cn (A/, A:) 


where the Crc and Vn functions are related to the Sn function by: 


Cn 2 x = 1 - Sn 2 x 
Vn 2 x — 1 — k 2 Sn 2 x 

The constants A and k are given by: 


A = H 


(Ji ~ [zM ~ 7 a) 


N 


{hhhD) 


(4.92) 


k = 


\ 


ih - w - g) 

(/1 - / 2 )(X> - / 2 ) 


0 < A- < 1 


(4.93) 


The above solution corresponds to a rotation about the / r axis with the initial 
condition u> 2 (0) = 0. 


For small values of k(I 2 ~ / 3 ) the elliptic functions approach the trigonometric 
functions: 


Sn(Xt, 0) = sin At 


Cn(\i, 0) = cos At 


2>n(At, 0) = 1 

The determination of the angular velocity components as functions of time does, of 
course, not complete the solution of the problem, because we still have to find the 
orientation of the body relative to an inertial observer at any time. This can be 
done, in principle, by substituting the solutions into Euler’s kinematical equations, 
which yields three differential equations for the three Euler angles. The solution of 
these equations in this general form presents a formidable problem and will not be 
pursued any further. However, it is important to notice, that although the angular 
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velocities are periodic functions of time, the motion of the rigid body as viewed 
from an inertial observer is no longer periodic. 

C) Perturbation Method 

In many practical applications the rigid body motion must not deviate greatly 
from a nominal direction. In such situations valuable insight into the dynamic 
behavior of a system can be obtained by linearizing the equations of motion about 
a nominal condition. 

To illustrate the method we consider an unsymmetrical rigid body which spins 
about the I r axis. If the asymmetry is rather small (i.e. / 2 » h) we observe from 
the Poinsot construction that the motion is confined to the vicinity of the spin axis. 
Therefore the first equation of Equation 4.82 is approximately 


IiLUi = 0 or u; x = uj 0 = CONSTANT 
The second and third of these equations can then be w’ritten as: 


cj 2 4" A^u ^3 — 0 

(J-9-1) 

U 3 — A 2 cj 2 — 0 

(-1.95) 


where 


Ai 



A 2 — 



u 0 


Because of the triangle inequality rule of the moment of inertia both | Ai | and | A 2 | 
are smaller than the nominal spin luq. 

We assume in the ensuing discussing that / 2 is the intermediate moment of inertia. 


In order to be able to compare the approximate solution with the exact solution 
obtained in the preceding section, we assume the same initial condition: 


u> 2 (0) = 0 for i = 0 

The equations will be solved by the Laplace transformation method. 


119 


Denoting the Laplace transforms of the angular velocities by capital letters: 


£M0} = «W 

we obtain: 


sfi2(®) “I - AifisW — 0 


(4.96) 


sfi 3 (s) — lj 3 (0) — A2fl2(s) — 0 (4.97) 

where u> 3 (0) is the initial angular velocity along the / 3 - axis. 

Solving for the Laplace transforms of the angular velocities we obtain 

The corresponding inverse transforms yield the angular velocities as a function of 
time: 


^2(0 = sin yJXiX^t (4.99) 

v Ai A 2 

u 3 (t) = ^ 3 ( 0 ) cos ^Aj A 2 < (4.100) 

For a spin close to the I x - axis the quantity D ^ J x . Inserting this approximation in 
Equation 4.93 it is seen that A = V% A 2 , that means the angular frequencies of the 
angular velocities are in agreement with the exact solution. It can also be easily 
verified that the same agreement holds true for the amplitude ratio of u; 3 and u; 2 . 

It is also of interest to plot the path of the angular velocity in the equatorial 
x? ~ 2*3 plane of the body. This is actually the projection of the polhode onto the 
equatorial plane. \Ye obtain an elliptical path which is traversing the arj-axis in a 
clockwise sense if > 0(7] > J 3 ) and in a counterclockwise sense if Aj < 0(/i < / 3 ). 
This is also in agreement with the Poinsot construction. 

The real interest lies, of course, in the motion of the body, especially the axis of 
symmetry or spin axis, as seen by an inertial observer. To this end we introduce the 
modern Euler angle system which is particularly suited for perturbation studies. 
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Repeating the corresponding kinematics] equations we have: 

u i = <p — V sin 8 

u) 2 = 8 cos <p + xp cos 8 sin <p 

u >3 = xp cos 8 cos (f> — 8 sin (p 

We now use the approximation u>i m = <p and restrict the motion to small angles 
8 and xp. Thus: 

uj 7 — 8 cos <jj 0 t + ip sin 
cj 3 = — 8 sin uj 0 t + xp cos 

CASE A: Spin about Minimum Moment of Inertia (Rod: \ x < 0 ; / L < / 2 < h) 

Set: 


A = 




C 3 = W 3 (0) 


The equations of motion are then: 


u> 2 = c 2 sin \t = 6 cos + $ sin iu 0 i 


u> 3 = c 3 cos \t = — 9 sin ui 0 t + t/> cosw 0 < (4.101) 

We now introduce a complex cone angle: 

a = ip + i 9 (4.102) 

Adding the two equations of motion in quadrature we obtain: 

w 3 4- iw 2 = c 3 cosA< + ic 2 sin \i = (ip + id)e ,uot = de ,uot (4.103) 

Solving for the complex cone angle yields: 

a — (c 3 cosAl + ic 2 sin A<)e _, " 0< (4.104) 

To bring this equation into a more convenient form for integration we introduce: 

c 3 = B 2 and c 2 = Bi — B 2 (4.105) 
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Thus: 


a = [(l?i -f B 2 ) cos Xt + i(Bi — B 2 ) sin Xt] e ,UJBt 

— [(fii(cos Xt + isin Xt) + £ 2 (cosA< - isinAf)] e~ ,w ° f (4.106) 
= B 1 e -,(u ' 0-A)< + £ 2 e _,(wo+A) ' 

Integrating with respect to time furnishes the complex cone angle: 


B i 

(w 0 - A) £ 


:(u/ 0 -A)f 


th 

(<^o + A) 


— t(wo + A)f 


Converting back to the original constants by observing from 4.105 that: 


(4.107) 


B i 


C3 + C 2 

2 


and fi 2 


c 3 


~ C2 

2 


we finally arrive at: 


q = i [.4 1 e -,(u ' 0-A) ' + .4 2 e~ ,(u ' 0+A) '] 

where 


(4.108) 


C3 + C 2 

2(o-o - A) 


and 


c 3 ~ C 2 

2(o) 0 + A) 


It is seen that the complex cone angle is represented by two rotating vectors having 
different amplitudes and frequencies. 

In order to make a more definite statement about the motion we have first to 
prove the following two inequalities: 


Proof# 1: 


A = yAiA^ < o-o and 



(4.109) 


A) A 2 


[h - h){h - ij) 

hh 


0-0 < 0-0 


{1\ — ~ I 2) < I2I3 
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I x — I \I 2 — I\I 3 + I2I3 < I2I3 

/i(/i -h-h) < 0 Q- E. D. 

The final inequality holds because of the triangle inequality for the moments of iner- 
tia. As a consequence both vectors in Equation 4.108 represent clockwise rotations. 

Proof #2: 

This inequality ensures us that I 2 is indeed the intermediate moment of inertia. 

h. — Hi _ lilh. > 1 
A 2 (A — h)h 

Because I x < I 3 and I x < I 2 the inequality can be rewritten as: 

(h ~ h)h > (h ~ h)h 

Adding the term / 2 / 3 to both sides of the inequality gives: 

[h + h ~ Ii)h > (h + h ~ h)1 2 

According to the triangle inequality rule I? + I 3 — I x > 0 and therefore, 

h > h Q- E - D - 

The motion can be made visible by projecting the tip of the unit vector along the 
body-fixed avaxis (spin axis) onto the inertial Y - Z plane which is normal to the 
nominal spin direction. 
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The coning motion is a superposition of a large amplitude small frequency rotation 
and a small amplitude high frequency rotation both having a clockwise rotation 
about the spin axis. The spin itself is, of course, also clockwise. The motion is 
confined within an annular ring. The minimum cone angle is the initial one: 


-T + 2 — &A//A' 


■^3^3(0) 

0 


The maximum value is obtained for (T 2 < 0!) 


-4, 


— A 2 — a.X 


hjh — 1 1 ) 



~ h) 


(*m 1 a 


This is in agreement with the exact solution on page 115 and 116. The motion is 
in general, not periodic. 


Case B: 


Spin about the Maximum Moment of Inertia (Disk: A ; >0:7] > / 2 > / 3 ) 

This case is obtained by changing the sign in the uj 7 equation of Equation 4.101 
because C 2 < 0 for A! > 0. All other steps remain the same. 

The final complex cone angle is: 

a = i [,4i e _,( " 0+A)t + .4 2 e _( “ 0_A)( j (4.110) 

where now 


-4, 


C 3 + C 2 , , C 3 - C 2 

— rr and A-> = — 

2(wo + A) 2(u 0 — A) 


The coning motion consists now of a superposition of large amplitude high frequency 
rotation and a small amplitude small frequency rotation both having again the same 
clockwise sense of rotation. This is graphically illustrated in the following figure 
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The motion is again confined to the annular region. Maximum and minimum cone 
angles are identical to the previous ones. 

NOTE: 

For a symmetric satellite / 2 = I 3 the coning motion becomes steady with a frequency 
of 


UJ = UJ 0 + \ = 7 - Wo Precession 

h 

This precessional motion is always clockwise (“forward ) about the spin axis 
regardless of whether the spin is about the maximum or minimum moment of in- 
ertia. One can obtain this result from the preceding steps by observing that for 
a symmetric body / 2 = h which leads to c 2 = c 3 . It is important to pay careful 
attention to the sign of A which changes when going from a long slender rod to the 
case of a flat disk. It is negative for the former case and positive for the latter. 

Sometimes the geometric axes or control axis system deviate slightly from the 
principal axis system. Let us assume that the control axis has a small angle 0 w ith 
the spin axis. This misalignment is strictly a geometric effect and does not enter the 
dynamics equations. It can be simply taken care of, by adding to the complex cone 
angle of Equation 4.108 or Equation 4.110 a vector of magnitude 0 which rotates 
clockwise with the spin frequency u) 0 , is an example . Equation 4.108 would read 
for this case , 

a = i [T 1 e~ ,( “' 0-A)< + T 2 e- ( “ 0+A) ' + 0e~' uat \ (4.111) 
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4.6 Forced Motion of a Rigid Body 


Unlike the case of the torque-free rigid body the equations of motion of a rigid body 
which is acted upon by external torques can only be solved analytically for very 
special cases. It was first shown by Lagrange (1788) that the rotational’ equations 
of motion can be integrated for a symmetric body in a uniform gravity field with a 
fixed point on its axis of symmetry (“heavy top”). We will not discuss this rather 
length} treatment but instead again will present perturbation methods which can 
be used to describe small deviations of a rigid body from a nominal condition. The 
mathematical models will be applicable for a variety of spinning (“gyroscopic”) 
systems. The methods will differ depending on whether or not the external torques 
depend on the orientation of the bodv. 

Method I (Body-fixed torques) 

In this case the Euler equations of motion can be directly integrated to find the 
angular velocities as functions of time. The motion relative to inertial space is then 
obtained by subsequently integrating the kinematical equations. The method is 
essentially the same as the one used in the preceding paragraph for the torque-free 
case. 

As an example we consider an axially symmetric spinning missile with a thrust 
misalignment. The principal axes system is aligned such that the thrust misalign- 
ment produces a constant torque about the £ 2 -axis. 

The equations of motion are then 


h u 1 = 0 


/? w 2 + (/, - J 3 ) un u> 3 = L 7 / 2 = / 3 


h W 3 + (h — I\)u>i w 2 = 0 

From the first of these equations we see that the angular velocity component w , 
about the spin axis is constant. We will denote it by u 0 . Introducing the quantity: 


h ~ h 

= — } w 0 
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we can write the other two equations as: 


CJ2 + A U3 — — 

h 

CJ3 — A U2 — 0 

This set of equations can be solved again by the method of Laplace transformation 
and yields the solution for the initial condition u; 2 (0) = 0: 


,j\ Au; 3 (0) . 

«a(0 = — rrr sin 


A | t + 
L, A 


I A | / 3 


sin 


w 3 (<) = w 3 (0) cos | A | i - cos I A | t + 


A | t 
\L 2 


where w 3 (0) is the initial angular velocity about the / 3 -axis. 

For the symmetric case, there is no need to distinguish mathematically between 
a spin about the minimum and a spin about the maximum moment of inertia 
because the change in sign of A can be automatically performed by the sine function. 
Therefore, we obtain for both cases: 


w 2 (0 = ~ |w 3 (0) - jy- sin \t 
w 3 (<) = [w 3 (0) - jj - ] cos At + 


L2 
A h 


Using the approximate kinematical equations on page 114 and the complex cone 
angle introduced in Equation 4.102 we get: 


d = L(0) - -f 1 + fy-e- 

XI 2J X12 

Upon integration the cone angle is obtained as 


— JWQ* 


. [^3(0) - - lW „ -1 — 

(u^o + A) A/2 UJq 


1 e -iu 0 t 


(4.112) 


(4.113) 


For a spin about the minimum moment of inertia (A < 0) the resultant coning 
motion consists of the regular low frequency precession with frequency ui p = (cj 0 + 
A) which is superimposed by a small amplitude high frequency oscillation. Both 
rotational vectors revolve about the spin axis in a clockwise sense. 
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Figure 4.4: Floating Coordinate System 


KQTE: 

This method can also be used for an unsymmetric rigid body as long as the 
asymmetry is relatively small (/ 2 = / 3 ). Otherwise, the deviation from the nominal 
spin condition do not remain small and the linearization process used in this method 
is no longer valid. 

METHOD II (Modified Euler Equations) 

For axially symmetric bodies a modification of Euler's equation is possible by 
introducing a reference frame which is aligned with the symmetry axis but does not 
rotate with the body (“floating 1 ' reference frame). Because of the symmetry this 
does not cause a time changing moment of inertia matrix. The floating coordinate 
system is defined by the first two rotation sequences of the classical Euler angle 
system. The corresponding rotation matrix and kinernatical equations are therefore 
obtained by setting <£> = 0. 

The floating coordinate system is also known as node-axis system. The nominal 
spin is about the / 3 -axis. 

We can perform the desired modification by going back to the vector-dyadic 
form of the Euler equations: 


I'W + wxI-u = L 0 (4.114) 

It is important to realize that the angular velocity appearing in this equation is that 
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of the rigid body relative to inertial space. This is now expressed as the sum: 

uj = w 0 + Wf (-4.1 15) 

where u> 0 is the spin rate of the rigid body and is the angular velocity of the 
floating frame. In vector form, we have: 


UJq — </> e 3 U?o — <f> £ 3 (^-1 16) 

u) F = 9 e t + ip sin 9 e? + ip cos 9 e 3 

Also for the moment of inertia dyadic we have: 

I = h e, e] + I 2 e 2 + I 3 c 3 e 3 with I\ = I 2 (d. 117) 

Introducing these terms in the Euler equation yields: 

I ■ (uj 0 + w 0 f) + (wo + Wf) x 2" • (cj 0 + u > F ) = L 0 (4.118) 

It is easily observed that the term: 


u>0 X I ■ u?o = 0 

Therefole the Equation 4.118 can also be written as: 

I ■ (u?o + tb/') + UJ/r x T ■ (u?o + Wf) = Lo (4.119) 

This is the desired modification of the Euler equations. They can be expressed in 
terms of the Euler angles by substituting Equation 4.116 and Equation 4.117 into 
Equation 4.119 and carrying out the vector-dyadic operations. 

Many cases of practical interest can be investigated using the modified Euler 
equations. Linearization is often possible about a steady state condition. When 
analyzing gyroscopic systems, it is usually assumed that the spin rate (f> is large in 
comparison to ip and 9. It is then approximately equal to u>o which is the constant 
nominal spin rate. 

Example: 

Consider the case of the heavy top where the external moment is given by 
L l = 171(7/ sin 0 (see figure). 
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Assuming a steady state condition of 9 = f we obtain the set of differential equa- 
tions: 


I\ ip — 1 3 wo 9 — 0 
7i 9 + I 3 u> o i> = rngl = L x 

The solution of this set is most conveniently obtained by the Laplace transform 
method. 

Introducing the notation 


P = 


^Wq 


the angles ip and 6 are then for V'(O) = 0 and 0(0) = 0: 


1 , • L i . . L x 

V = - (V»o - ~j ) sin pi + i 

P -*3^0 * 3 ^0 


9 = - (wo - cos p t + ^ - - (i' 0 - 

P /3W0 2 p /3W0 

The gyroscope rotates about the z-axis with average angular precession 


u . p = 


u 

I3 Wo 


METHOD III: (Direction-dependent torque) 

In many practical cases the external torques acting on a spinning body depend 
only on its direction relative to a nominal spin direction. The equations of motion 
are derived in terms of the deviation of the rigid body from the nominal steady 
state spin condition. The body-fixed reference frame is aligned with the principal 
axes and the nominal spin rate w 0 is aligned with the inertially fixed x-axis having 
unit vector i. The total angular velocity w of the rigid body is again considered to 
be the sum of two components: 


W — Wo + U>B 
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where u) 0 is the nominal spin rate and u B the angular velocity of the rigid body 
associated with the deviation from the nominal reference frame which itself rotates 
with constant angular velocity ut 0 relative to inertial space. 

The nominal spin vector is expressed in terms of the deviation angles of the 
rigid body from the nominal reference frame. Using modern Euler angles we obtain, 
therefore: 


uj 0 = aj 0 i = w 0 [(cos tp cos 0)ei + (— sin rp cos <p + cos ip sin 9 sin <p) e 2 

+(sin ip sin <p + cos ip sin 9 cos <^>)e 3 ] 

The angular velocity of the body frame relative to the nominal frame is given by 
the kinematical equations: 

uj b = - xp sin 0)e! + (0 cos <p + ip cos 9 sin <p)e 2 

+(ip cos 9 cos <p — 9 sin <^)e 3 

As mentioned before, the reference frame is aligned with the principal axis sys- 
tem which means the moment of inertia dyadic is given by 

X = I\ ei ei + J 2 62 62 + 7 3 e 3 e 3 (4.120) 

Euler’s equations are then: 

I • (w 0 + u B ) + (u>0 + u B ) X I ■ (w 0 + U ) B ) = L (-1-121) 

Carrying out the various vector-dyadic operations and retaining only linear terms 
results in the following set of perturbation equations. 


h4> = L x 


I 2 9 + (/1 — / 3 ) ul 9 — (I 2 + 1$ — I\) uq \p — Z-2 


I 3 xp + ( I\ — h) Wq ip + {h + h ~ h) w 0 1 9 — 7-3 

It was assumed that the deviations of the rigid body from the nominal spin state 
remain small. We can see from the first equation, that the roll angle <p will increase 
with time if there is a torque about the 7i-axis. For the equations to be valid we 
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have to assume, therefore, that L x = 0. Furthermore, for an unsymmetrical body, 
the deviations of the spin axis from the nominal direction are only small if the 
asymmetry is small i.e. , / 2 % J 3 . 


The following examples will be restricted to the symmetrical condition / 2 = / 3 . 
NOTE: 

This method can also be used in the absence of external torque (torque-free 
rigid body). In this case, Method 1 is superior because it works with two pairs of 
first order differential equations rather than with a pair of second-order differential 
equations. The results are, of course, identical. 

Example 1: Heavy Top. (/ 2 = / 3 ) 

e a 

COo 


Introducing the notation: 

K = (/, - / 2 ) L ol - mgl 

G = 2 h - I, > 0 

we can write the perturbation equations as: 

/ 2 9 + K 9 — G ljq tI> = () 

1 2 ^ T K xjj + G u>q 9 = 0 
Introducing also the complex deviation angle 

8 — xp T i9 

we can add the Equation 4.122 in quadrature to obtain 

1 7 8 - i G u> 0 8 + K 8 = 0 (4.124) 
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(4.122) 

(4.123) 




The external torque is: 

L 0 = 1 x W 

= / e, x (— mgi) 
i = ei - ti> e 2 + 9 e 3 
L 0 = mgl 9 e 2 4- mgl xp e 3 
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Solving by Laplace transformation we have the characteristic equation 


The characteristic roots (eigenvalues) are: 

Gu! 0 


> 1,2 — 1 


2 /, 


G ujq s k — 0 


(4.125) 

are: 



:± \Z r 2 p !+ 

K 

h. 

(4.126) 


or in terms of the original system parameters: 


s 1,2 — 1 


(2/ 2 -/i)w 0 , fj lW0v 2 rng[ 

2h V 2 1 7 ] h 


— i-^1,2 


(4.127) 


The complex deviation angle can perform a uniform rotation with frequency A! or 
A 2 (or both) as seen by an observer moving with the rigid body. Rotational motion 
exists, however, only as long as the radicand is positive. This means that for a 
stable motion, the spin rate has to be higher than 


ljo > j- \J rngl / 2 (4.128) 

h 

In this case, we have the so called “sleeping” top. Of greater interest is the rotational 
speed of the deviation vector as seen by an inertial observer. To do this we just 
have to add the clockwise rotation of the nominal reference frame. The complex 
deviation angle in the inertial frame is therefore: 


6j = 6 0 e i(A - wo)< (4.1 29) 

where 6 0 is the initial condition. We obtain then two “precessional” frequencies, a 
low one and a high one. 


u>x = - 



LOW 


(4.130) 


i w o //7 i^Qs? 

2/7 "y 2/ 2 h 


HIGH 


(4.131) 


According to the definition of the complex deviation angle given in Equation 4.122 it 
can be observed that both precessional frequencies indicate a clockwise ( forward ) 
precession. 
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The precession of the heavy top is usually the slow precession. It is interesting 
to notice, that the high precession rate emanates from the precession rate of the 
torque-free symmetric body as given on page 119. 

If the spin rate is high we can approximate the low root by expanding the square 
root in a Taylor series and obtain 


or in another form: 


<*>i 


mgl 
h w o 


(4.132) 


mgl — ui\ I ujq 


(4.133) 


Example 2: Spin-Stabilization 

Artillery shells (bullets) and missiles which are exposed to destabilizing aero- 
dynamic forces are often stabilized by giving them a high spin rate about the long 
axis. 



In this case the destabilizing moment is due to the resultant aerodynamic force R 
acting at the center of pressure (C.P.) which is in front of the mass center (C.M.). 

The torques exerted by this force about the transverse axes can be derived in a 
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similar manner as in example 1 and yield: 


L = a j£ )te 62 + {q3 ~^) t + 63 

where q = j p v 2 is the dynamic pressure and s a reference area, and the 

normal force coefficient. 

It is seen that this example is identical with the heavy top. All we have to do is to 
replace 


dC N . 

mgt - 


Example 3: Rotor Dynamics. 

An important aspect of rotor dynamics is the gyroscopic motion in which the rotor 
system is precessing in its bearings. Assuming a symmetric rotor and bearing 
support, the elastic and dynamic characteristics depend only on the direction of 
the spin axis relative to its nominal state. The effect of the bearing stiffness is 
considered to be equivalent to a linear spring, with spring constant k. Unlike the 
previous two examples the spring forces are restoring forces. 


z 



Figure 4.5: Rotor Dynamics 
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The equations of motion are again identical with Equation 4.122 if we set 


k = (I 1 - I 7 ) u 7 + 2 k e 


In rotor dynamics jargon the precessional motion of the rotor is referred to as 
“subsynchronous whirl.” The two whirl frequencies are then: 


w 1 = - 


Id 2 — — 


. 2 kP 

irw + ~ir 

1^0 , 2kP 

wT + V'w’ + "77 


LOW 

HIGH 


(4.134) 

(4.135) 


The low frequency precession is always counterclockwise (“backward whirl”) whereas 
the high frequency precession is clockwise or forward. Since the inevitable unbal- 
ances in a rotor have the tendency to excite a forward precession, the backward 
whirl is usually not observed. For the heavy top both precessions are forward and 
only the low frequency is usually observed because it requires less energy for it being 
excited. 


NOTE: 

The present model allows a quick and rigorous assessment of the stability behav- 
ior of a gyroscopic system in the presence of damping by invoking the Kelvin-Tait 
Theorem (1921) . According to it internal damping destroys gyroscopic stability if 
the stiffness matrix of the linearized differential equations written in matrix form 
is negative definite. Therefore, a torque-free satellite spinning about the minimum 
moment of inertia is unstable when damping is present. For the rotor dynamics 
example inspection of the stiffness matrix shows that the system becomes unstable 
if the rotor speed is 


2 2 k l 7 

u o> I — j 

This phenomenon of instability due to internal friction (damping) is well attested. 
Method IV: Approximate Theory 

If the spin rate u> 0 of a symmetric rigid body is very high it is possible to derive a 
simplified model of its gyroscopic behavior which is sufficiently accurate to explain 
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many practical applications. Assuming that the spin w 0 is practically constant the 
Euler torque I • w in Euler’s equation can be neglected. We obtain then: 

L 0 = w x H where H = ! ■ w (4.136) 

We also assume that the angular momentum is essentially along the spin axis, i.e. 

H = I 1 u? 0 (4.137) 

Since the magnitude of the angular momentum is constant, it can only change its 
direction which, of course, is now also identically with the directional change of the 
spin axes. This directional change is called the “precession” of the gyroscope and 
is governed by Equation 4.136, where the angular velocity becomes the precession 
rate u> p . We also still observe that the angular momentum change is equal to the 
applied torque if the reference point is fixed or at the mass center of the gyroscope. 
Equation 4.136 is then written in the form: 

L 0 = ^ = w p x H = h (w p x w 0 ) (4.138) 

at 

Another formulation which aids the understanding of the gyroscopic behavior is: 

dU = L 0 dt (4.139) 

This means that the change of angular momentum vector is in the direction of the 
moment applied. 

A few examples are presented to illustrate the applications of these considera- 
tions. 

Example 1 Bicycle wheel. 

Consider the front wheel of a bicycle. The angular moment of the wheel is to the 
left. 
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If the rider tilts to the right, the external moment L is positive along the forward 
speed direction and causes a precession of the wheel to the right. The ensuing 
centrifugal force restores the bicycle to its upright position allowing a free-hand 
ride. 



The friction on the surface produces a positive torque about the C.M. in the upward 
direction. It has a tendency to bring the heavy top in the upright position. 

Example 3 Tippy Top (Class Ring) 


/ H 





The friction force being directed out of the paper plane produces a downward torque 
which precesses the top such as to invert its orientation. 
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4.7 Rheonomic Systems 


Dynamic systems containing internal mass elements whose motions relative to the 
main body are explicit time functions are known as rheonomic systems. It is appar- 
ent when inspecting the translational Equation 3.4 and the rotational Equation 3.12 
of section III that such internal moving parts affect the overall dynamics of the main 
body. If these effects are undesirable they induce disturbance forces and torques 
(e.g. crew motion, running machinery, etc.). If these motions are used to produce 
desirable effects they result in control forces and torques. In this case, the mo- 
tion is usually restricted to rotational motion. We will therefore refer to those as 
gyroscopic systems. 

For simplicity, the following discussion will deal only with a single internal gy- 
roscopic system. Generalization to a situation involving multiple internal moving 
elements should be straightforward. 



Figure 4.6: Internal Gyroscopic System 

Let fi be the angular velocity of the main body Mo relative to inertial space and 
be the angular velocity of the rotating mass m relative to the main body. If the 
origin of the main body reference frame coincides with the total mass center then 
the rotational equations are: 




i 


R x [2 (/? x v) + a ] dm = L 


(4.140) 
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where 1 is the moment of inertia of the total system including the moving mass. In 
most applications it can be assumed that the change in orientation of the rotating 
mass caused by its relative motion has a negligible effect on the overall moment of 
inertia, which is therefore, usually treated as a constant. This must, of course, not 
necessarily be true of all systems. 

For the rotational motion of the internal mass, we have the relative velocity and 
acceleration: 

v = « xr 

a = wxr + wx(wxr) (4.141) 

where r is the position of the moving mass with respect to its own center of mass. 
Therefore: 

J r dm = 0 (4.142) 

The position vector R of a moving mass element can be written as: 

R = ^ + r (4.143) 

where & is the position of the mass center of the moving mass relative to the mass 
center of the total system. 

We can write then the Coriolis term in Equation 4.140 in the form: 

2 J R x (/? x v) dm = 2 J (£ + r) x [fl x (u> x r)j dm 
— 2 J r x [f2 x (u> x r)] d m = 2 j (f2 • r)(r x u>) dm 
Introducing the Coriolis dyadic 

C m = J (r r) d m = (^ tr I m )£ - I m (4.144) 

where X m is the moment of inertia dyadic of the moving mass. The Coriolis term is 
then: 

2 J R X (n X v) d m = 2/2 • C m x u> (4.145) 

The relative acceleration term in Equation 4.140 can also be written in a more 
convenient form by using Equations 4.141 and 4.142: 
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(4.146) 


J R x a d m = |rx[wx(«xr) + (wxr)] dm 
= J m ■ u > 4 w x I m • U7 

The final form of the rotational equation for the rheonomic systems is: 

l-n + f2xl-n + 2f2 C m xw 

+ I m • u> + u>xI m -u> = L (4.147) 

a) Gyroscopic Attitude Control 

There exist in practice essentially two different techniques of controlling the 
orientation of a spacecraft in space using internal gyroscopic devices. These devices 
are generically called (angular) momentum exchange controllers . They are both 
using the two relative acceleration terms of Equation 4.147 to affect the desired 
attitude control. For the purpose of attitude control analysis and simulation, these 
terms are brought over to the right hand side of Equation 4.147 and then looked 
upon as effective control torques about the mass center of the main body. 

1) Momentum Wheel 

This device is also known as reaction wheel, inertia wheel or flywheel. In this 
technique, the spin axis of the controller is kept in a fixed direction relative to the 
main body and only the magnitude of the rotational speed of the spinning \\ heel is 
changed. Since the wheel is spinning about a principal axis the term u) x T m ■ u> is 
zero. The so called control torque of the device is then: 

L„ = -Z m - u> c (4.148) 

where u> c is the controlled change of the spin rate which is determined bj the desired 
control torque L c . 

2) Control Moment Gyro (C. M. G.) 

In this technique the rotational speed of the device is kept constant at a very high 
level. The desired control torque is generated by a proper change of the direction 
of the spin axis. This change can be affected by a single and two degree of freedom 
(single and double gimbal) C. M. G. system. 

Since the spin rate of the C. M. G. is very high and practically constant, the 
term 7 m • u> in Equation 4.147 can be safely neglected. Furthermore, the angular 
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momentum vector can be assumed to coincide with the spin axis of the C. M. G. 
rotor. Therefore: 


H m = I m ■ u> = I ro • u?o (4.149) 

where u> 0 = spin rate vector. 

The control torque of the C. M. G. can then be approximated with a very high 
degree of precision by: 


L c = u? c x H m (4.150) 

This particular form of the C. M. G. control torque allows a highly useful treatment 
of the C. M. G. behavior in terms of angular momentum considerations especially 
when, as in most practical applications, multiple C. M. G. systems are employed 
for attitude control. 


A major task in the proper use of C. M. G.’s is to find optimum “steering” laws 
for the directional changes of the C. M. G. spin axes and to avoid “saturation” of 
the C. M. G.’s, a condition in which the controllers loose their control effectiveness 
along certain directions. Some of the difficulties associated with this task stem 
from the fact that the solution of Equation 4.150 in terms of the unknown angular 
velocities u? c is not unique, because it involves, in general, more unknowns than 
equations. 

The literature in this area known as angular momentum management is quite 
extensive. Any further discussion would go beyond the scope of the present expo- 
sition. 

m m 

No viable implementation is known which combines both techniques in one con- 
troller. 

b) Gyroscopic Reaction Torques 

In many practical applications it is required to determine the reaction torques 
which a gyroscopic system exerts on its surroundings (e. g. bearings) when its 
motion is known. The solution of this problem is relatively simple. Knowing the 
motion of the gyroscope one can calculate the left-hand side of Euler’s equation to 
obtain the external torque which must be applied to generate the assumed motion. 
The reaction torque is then represented by an equal torque of opposite sign. 
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If there exists, however, a considerable interaction effect between the gyroscopic 
system and the main body it would be necessary to work with the rotational equa- 
tion (4.146). Application of this equation requires the evaluation of the Coriolis 
term. Transferring this term to the right-hand side yields the reaction torque aris- 
ing from it as: 


L r = -2J7 -C m xu> (4.151) 

If the main body angular velocity is resolved along the principal body axis system 
of the gyroscope corresponding to: 


Q — + Q 2 e 2 T 

the Coriolis reaction torque becomes in component form: 

Li = Q* w 2 (h + h ~ h)~~ ^2 w 3 (fi — ^2 + h) 

L 2 = fil U >3 ( 1 2 + I 3 ~ h)~ ^3 w i (^1 + I? ~ h) 

L3 = Q 2 <*>i (h ~ h + h)~ u 2 {h + — h) 

Reaction torques can, of course, also arise from the other relative motion terms. 
Example 1: Airplane propeller. 

Determine the gyroscopic reaction torque of a two-blade airplane propeller when 
an airplane makes a horizontal turn with constant yaw rate ip. Rotational speed of 
the propeller is u> 0 . 



Figure 4.7: Airplane propeller 
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The moment of inertia of the propeller is approximated by a uniform slender rod of 
length i: 


I = h = h 


m l 2 

~W 


The pertinent angular velocities are: 


h = 0 


Q 2 = \i) sin <p CI 3 = ipcos <p Qi — 0 

U>\ LJq U>2 — 3 — 0 

The Coriolis reaction torque becomes then: 


Ti=0 X 2 = —2 1 ip luq cos <p L 3 = 0 

This torque can be referred to the airplane coordinate system (X, V, Z) as: 

Ly = — 2 / ip u>o cos 2 (f> 

L% = —2 I xp lj 0 sin <p cos (p 

where (p — / 

Notice that the frequency of this torque is twice the frequency of the propeller 
angular speed. By the way a three-blade propeller (X 2 = J 3 ) eliminates this cyclic 
propeller torque. Since the propeller is spinning about a principal moment of inertia 
at a constant speed no reaction torques arise from the last two relative motion terms 
of Equation 4.1-47 However, there is an interesting torque coming from the second 
(centrifugal) torque term. Remember that the moment of inertia in this term refers 
to the total moment of inertia of the main body and the moving mass. Looking 
only at the torque coming from the propeller we have (J p = propeller moment of 
inertia). 


L fi - —fi x l p • fi 

= +fl 2 f 1-2 e t = / 2 i/> 2 sin <p cos <p 

This torque is along the x r axis which coincides with the inertial x-axis. The cor- 
responding component is: 


L x = / 2 ip 2 sin <p cos <p 
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It is not a gyroscopic reaction torque because it really depends only on the orien- 
tation of the propeller relative to the airplane and exists also when the propeller 
is not spinning. In fact, it is related to the internal shear stresses which are even 
present in a rotating system. The torque disappears for (f> = 0 and <f> = 90° when the 
airplane rotates the propeller about a principal axis and therefore does not induce 
a centrifugal torque. 

NOTE: 

Euler’s equations could be used directly to determine the reaction torques. The 
present method shows more clearly how an internal gyroscope affects the total 
systems dynamics. It is seen that the gyroscopic reaction torque of a symmetric gy- 
roscopic system is exactly given by the Coriolis torque of Equation 4.151. Assuming 
a steady spin about the I \- axis such that: 


uj = u> 0 = u> 0 ei u> 0 = CONSTANT 

and axial symmetry, which means / 2 = h, we obtain the gyroscopic reaction torque 
of a symmetric gyroscope from the components of Equation 4.151 which reads then 
in vector form: 


L/i = — O3 ujq I ie 2 + 0 2 wo /l ©3 


or in another form: 


Lfi = -fix(/ 1 u 0 ) = -flxH (4.152) 

where H = I\ u>o is the angular momentum of the gyroscope relative to the rotat- 
ing total system. The negative sign indicates that the gyroscopic reaction torque 
always tends to bring the two vectors fl and u> 0 into coincidence. This is the 
principle of homologous parallelism as enunciated by Leon Foucault (1819-1868). 

It is very important to notice that the relationship established in Equation 4.152 
is an exact one and holds true regardless of whether Wo is small or large relative to fl. 
This is in distinct contrast to the approximate relationship given in Equation 4.138 
for the precessional motion of a symmetric gyroscope under the influence of an 
externally applied torque. Here we have the situation of an externally applied 
(forced) precession. 


145 


Example 2: Bicycle Wheel 


Consider the top view of the front wheel of a bicycle 



If the rider forces a precession on the front wheel by turning the handlebar to the 
left, the corresponding reaction torque will tilt the bicycle to the right according to 
the principle of homologous parallelism. (Foucault’s Principle). This maneuver can 
be used to initiate a right-hand turn. 

Example 3: Ore Crusher 



The angular velocity Q about the vertical axis causes each roller to have a relative 
angular velocity: 


Q £ 


LUq — 


r 
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The gyroscopic torque is obtained from Equation 4.152 as 


I L r |= h 


0 2 £ 

r 


By Foucault’s principle it will increase the normal force between the rollers and the 
surface. 
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Chapter 5 

Lagrangian Dynamics 


Newtonian or vectorial mechanics bases everything on two fundamental vectors: 
force and acceleration. Each component of a dynamical system is isolated and 
treated individually. Free body diagrams are introduced by which all forces are 
represented which contribute to the acceleration of the individual components. 
These forces include both the externally applied forces (known) and the internal 
reaction forces (unknown) acting on the isolated component. The latter are also 
called constraint forces because they act as to maintain the kinematical constraint 
conditions existing between the individual components of a dynamical system. To 
distinguish between the two forces one can call the applied forces also forces of 
physical origin” and the reaction forces also “forces of geometric origin” because 
they emerge from the geometric configuration of the system. Applied forces reveal 
their origin by the fact that their mathematical expression contain quantities which 
can only be determined by an experiment. 

Example: 

Static friction is a constraint (reaction) force whereas sliding (kinetic) friction is an 
applied force involving the experimentally determined coefficient of friction p. 

Lagrangian or analytical mechanics bases everything on tw'o fundamental scalars, 
kinetic energy and work. From a philosophical standpoint it is indeed surprising 
that two scalar quantities contain all the information regarding the motion of the 
most complicated system because motion is by its very nature a vectorial quantity. 


A component is no longer considered an isolated unit but a part of an overall 
system of interacting components. The great superiority of the Lagrangian dynam- 
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ICS over the vectorial dynamics stems mainly from the fact that knowledge of the 
unknown constraint forces is not required, but only knowledge of the kinematical 
conditions. Working with scalars rather than with vectors is another contributing 
factor to the simplicity of the analytical treatment. 

Before deriving the Lagrangian equations of motion it is necessary to introduce 
and define some distinctive concepts of analytical mechanics. 


5.1 Constraint Equations 

The kinematical constraints existing between systems components can be mathe- 
matically expressed as equations connecting conditions between their positions or 
between their velocities. 

Example: 


) 


Figure 5.1: Simple Pendulum 

Consider the simple pendulum constrained to move in a vertical plane. 
The position of the mass m can be defined as: 

r = xi + yj 

and its velocity as: 

v = xi + yj 


(5.1) 

(5.2) 
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One way of expressing the kinematical constraint is to state that the velocity of the 
mass particle must always be perpendicular to its position vector i.e. 

v • r = 0 (5.3) 

In component form Equation 5.3 reads: 

x x + y y = 0 (5.4) 

This is a kinematical condition on the velocity of the mass particle. It can be 
integrated and yields: 


X 7 + y 2 = e (5.5) 

This is a condition on the position of the mass particle and could have been obtained 
directly by observing that the pendulum mass is constrained to move on a circle. 

NOTE: 

It is often easier to determine the kinematical condition on the velocity of system 
components than on their position. 

If the conditions on the velocities can be integrated to yield conditions on the posi- 
tions of the system components the system (or the constraint) is called holonomic 
(“holonomic” is the Greek word for the Latin word “integrable”). If this integration 
cannot be performed the system (or the constraint) is called nonhqlQnQIIlic. 

Example: 

A characteristic and often quoted example is that of a ball which rolls without 
slipping on a horizontal plane. The kinematical constraint of “rolling (and piv- 
oting) requires that the instantaneous axis of rotation goes through the point of 
contact 0. 


U ) = + UI 2 C 2 + ^ 1^1 

v = + v 2 e 2 -I- v 3 e 3 

Ro = Ro^2 


151 



e 2 



Figure 5.2: Rolling Ball 

The velocity of the center C of the ball is then given by 

V = U X R 0 R 0 = R 0 e 2 (5.6) 

where R 0 defines the position of C relative to the contact point 0. Using a body- 
fixed reference frame with unit vectors e!,e 2 and e 3 , Equation 5.6 can be written 
in component form (“ROLL” CONDITION): 

= Ro l>2 = 0 l>2 = LJ\ Rq (5.7) 

These are three constraint conditions on the velocity of the center of the ball, but 
only one, namely v 2 = 0, can be expressed as a condition on the position. In 
fact it states simply that the ball has to move in a horizontal plane such that its 
center has a constant distance from the surface. The other two conditions cannot be 
integrated without solving the entire dynamical problem. This is due to the fact that 
the angular velocity u> of the ball is expressed in a moving reference frame whose 
orientation is given by Euler s kinematical differential equations of the preceding 
section I\ and this instantaneous orientation changes with time. It is, however, 
possible to establish differential equations of constraint from Equation 5.7 either 
in terms of Euler angles or by considering the angular velocity components , u> 2 
and u) 3 to be time derivatives of angular displacements a 1) a 2 , and a 3 (so-called 
quasicoordinates) about the instantaneous body axes such that: 
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^3 ~ 


(5.8) 


With these defined, 


w 1 


dai 

dt 



da 3 


the roll condition 5.7 can be written as: 


dx 1 + R 0 da 3 = 0 ; ( 5 . 9 ) 

d x 3 + R 0 d a x = 0 (5.10) 

But these differential relations are not integrable. 

A system for which the kinematical conditions change with time is called rheonomir 
otherwise it is called scleronomir . 

Examples of a rheonomic system is a mass particle which moves on a surface 
which itself is moving according to a prescribed time function. Another example is 
a pendulum whose length is a given time-function. 

The essential difference between rheonomic and scleronomic constraints is that 
rheonomic constraints do work. As a consequence rheonomic systems are not con- 
servative. This is the reason why rheonomic systems can become unstable in a very 
unsuspecting way. 

Illustration: 


Consider a tennis racket. If it is held fixed, the ball is reflected without change 
in energy. If the racket yields energy is taken out of the ball, if it moves against the 
ball, it transfers energy to the ball. 

Generalized Coordinates 

Lagrangian dynamics extensively uses coordinates other than Cartesian coordi- 
nates, which are then called generalized coordinates. They are any set of parameters 
which can be used to define the configuration of a dynamical system. Some of the 
generalized coordinates may not have geometrical significance and are therefore also 
called hybrid coordinates. For example, the amplitudes in a Fourier expansion of 
the position vector R. may be used as generalized coordinates. 

Using generalized coordinates ?i, < 72 , ■ • • , q n the holonomic kinematical constraints 
can be mathematically expressed as: 


=0 (j = 1 , 


2 , ...,m) 


(5.11) 
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Nonholonomic constraints are conventionally written in the form: 


£ Ci,(,ul 2 q„,t)d qk + c,(t)=0 (J = 1,2 m) (5.12) 

fc=l 

Degrees of Freedom (D. O. F.) 

The number of degrees of freedom of a system is equal to the number of independent 
generalized coordinates necessary to define the configuration of a dynamical system. 
This number is characteristic of a given dynamical system. If N parameters are nec- 
essary and sufficient to define the system configuration we sa^ it has N degrees of 
freedom.” Each independent kinematical constraint condition reduces the number 
of degrees of freedom by one. 


Examples: 

One D. O. F.: 
Two D. 0. F.: 
Three D. 0. F. 
Five D. O. F.: 
Six D. O. F.: 


A mass particle moving along a given curve. 

A mass particle moving on a given surface. 

A mass particle moving freely in space. 

Two mass particles connected by a massless rod (Dumbbell). 
A rigid body moving freely in space. 


NOTE 1: 

A ball (coin) rolling on a horizontal plane is sometimes said to have five finite 
D. 0. F. ’s and three infinitesimal D. 0. F. ’s. 


NOTE 2: 

Linear and angular velocity components of a rotating reference frame cannot be in- 
tegrated to furnish position and orientation. They are called nonholonomic velocity 
components. 

NOTE 3: 

Some authors only refer to independent coordinates as generalized coordinates. In 
this case, the number of generalized coordinates is equal to the number of D. 0. F.'s. 
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Virtual Displacement/Virtual Work 


A virtual displacement is an infinitesimal change of a generalized coordinate 
which is compatible with the kinematical constraints existing at that instant of 
time. Any moving constraints are temporarily stopped. This process is a kind of 
mathematical thought experiment. To emphasize its virtual character Lagrange 
introduced the special symbol 6. It has the usual properties of the ordinary differ- 
ential d; for example 6(sin 9) = cos 9 6 9. A virtual displacement can, of course, 
be applied to several coordinates simultaneously, but they are, in general, not in- 
dependent because of the kinematical constraints. 


The work done by a force F during a virtual displacement 6r is called virtual work 
defined by the scalar product: 


6 ir = F-6r (5.13) 

To illustrate the difference between actual work and virtual work consider the work 
done by the Coriolis force acting on a mass particle. 

d\Y = - 2m(u» x v)-rfr (5. 14) 

where dv is an actual displacement and is therefore related to the velocity as dr = 
vdt. Equation 5.14 becomes then: 

dW = -2m(w x v) ■ vdt = 0 (5.15) 

The actual work of the Coriolis force during an actual infinitesimal displacement is 
zero. The virtual work, however, is: 

6JV = -2m(u) X v) • 6r ^ 0 (5.16) 
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5.2 Principle of Virtual Work (Bernoulli 1717) 


Consider a mechanical system in motion. According to d’Alembert’s Principle, the 
vector sum of all forces acting on each particle is zero. These forces include the ex- 
ternal forces F;, the inertial forces I; due to the acceleration of each particle and the 
constraint (reaction) forces R, which maintains the given kinematical constraints. 
Thus, the equilibrium of the force system on each particle is: 

F, 4- 1, + R, = 0 where I, = — m,R, (5-17) 

The virtual work of all these forces during the virtual displacements must likewise 
be zero: 


SW = E(F. + 1, + Hi) • = 0 (5.18) 

The principle of virtual work establishes now the following postulate : 

“For any mechanical system, the virtual work of the constraint forces is zero.” 
Mathematically expressed it is: 

mt = Y, R . ^ r - = 0 ( 5 - 19 ) 

This postulate is the corner stone of analytical mechanics. Many scientists Ref. 
(2,3) consider it to be an additional axiom of mechanics which cannot be derived 
from Newton’s laws. Adopting this viewpoint, analytical mechanics is more than 
just a different mathematical formulation of Newtonian mechanics. 


5.3 Generalized Forces 


Another fundamental concept of analytical dynamics is that of the generalized force . 
The generalized forces acting on a dynamical system are determined b^ calculating 
the virtual work done by the external forces during the virtual displacements 6q t of 
the coordinates q,. Each virtual displacement 6</ t will produce the virtual work 

MV, = Q, % (5-20) 

where Qi is a quantity containing the external forces acting on the system. This 
quantity is called the generalized force Q, associated with the generalized coordinate 

Qi- 
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Figure 5.3: Double Pendulum 

To illustrate the situation, we consider a double pendulum which can be defined 
by the two independent generalized coordinates 8 X and 0 2 . 

To calculate the virtual work done by the external gravity force, we first define 
the position of the two mass particles m x and m 2 in Cartesian coordinates as follows: 

mass mi : Xi = £ x sin 9 x 

Vi = cos 8 X 

mass m 2 : z 2 = £ x sin 8 X + £ 2 sin 9 2 
y 2 = £ x cos 9 X + £ 2 cos 8 2 

Next we perform the virtual displacements of the two masses: 

6 x x = £ x cos 9 X 6 9 X 

6 y x = — £ x sin 9 X 6 9 X 

6 x 2 = £ x cos 8 X 6 8 X + £2 cos d 2 & 8 2 

6 y 2 = —£ x sin 8 X 6 8 X — £ 2 sin d 2 6 0 2 
The virtual work associated wdth these virtual displacements is: 


Ml — F x ! <5 :ei + F x2 6 x 2 + F yX Sy x + Fy 2 8y 2 
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Because the external forces are only due to gravity we have 

F x i - F x 2 = 0 and F yl = m x g , F y2 = m 2 g 
The virtual work is then: 

6 IT = m i g(—£i sin 6 X 8 9\) + m 2 g(—£i sin 8 X 88 x — £ 2 sin 8 2 8 6?) 

Collecting terms with 8 8 X and 8 8 2 yields: 

<5VT = —(mi + m 2 )g £ x sin 8 X 8 8 1 
—m 2 gl 2 sin 8 2 8 8 2 

According to Equation 5.20 the generalized forces associated with the generalized 
coordinates are then: 

for 8 ! : Q i = -(mi -I- m 2 )g sin 8i 

for 8 2 : Q 2 = —m 2 g £ 2 sin 8 2 

One could, of course, derive these two generalized forces more directly by finding 
the work done by the gravity force on the masses during an independent virtual 
displacement of their corresponding generalized coordinates 8 X and 8 2 . 

For the general definition of the generalized forces, we consider a system of N 
particles whose positions are defined by the vector r, and acted upon by the forces 

F,- 

The virtual work is then: 

r= £ Ff ■ i = 1 , 2 . . . A r (5.21) 

The positions of the particles are related to the generalized coordinates by: 

r . = Ti(qi,q 7 ...q m ,i) (5.22) 

The virtual displacements of the generalized coordinates are obtained by differen- 
tiating Equation 5.22 keeping in mind that a virtual displacement requires that 
8t = 0: 
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6r '~^ dq k 6qk ' * = 1 ’ 2 "-’ n 
Substituting Equation 5.23 into Equation 5.21 yields: 

* = 1 k = 1 

Changing the order of summation yields: 

A = 1 1 = 1 i 

The generalized forces are therefore: 


Qk ~ S F< ' Wk * ==1,2 ’"-’ n (5-26) 

The forces F, can be separated into external forces and constraint forces: 

F ' = F S £) + R. (5.27) 

Consequently, the generalized forces can be separated into generalized external 
lorces and generalized constraint (reaction) forces: 


= if< 


(£) d r, 


(5.28) 


Ql R> = 


(5.29) 


5.4 Classical Lagrange Equations 

Consider a mechanical system of N particles of constant mass m, with position 
vectors r,- in an inertial reference frame. 

The kinetic energy is: 


r= 2 ^ * = 1,2,..., TV (5.30) 

The kinetic energy can be expressed in terms of generalized coordinates using the 
transformation equations: 
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(5.31) 


r,- — r,(< 7 i , q?, • • • i <7n> 0 

Differentiating Equation 5.31 with respect to time gives: 


" dr, . dr , 
r ‘ = £ «■ 


(5.32) 


The kinetic energy is by substituting Equation 5.32 in Equation 5.30: 


N 


1 * A dr i . dr, 2 

r = 5 ? m ' ( f, dT„ , ‘ + ad 

* t=i fc=l ** 


(5.33) 


Define the generalized momentum p k associated with the generalized coordinate q k : 

dT . dTj t~ , 4 x 

p k = — = 2^ m,r, • — t 0 " 5 ^ 

dg* fr[ dq k 

The time rate of change of p k is: 


dp k ddT 


dr, 


dr. 


dt dt dq k ! 

Next using Equation 5.30 we calculate 

d T 


= - 37 (— ) = X>> r - • flT + m > r ’ • diJ 


dq k 

dr, 


— = Y. m, r, • 

dq k oq k 

This equation is seen to be identical with the last term of Equation 5.35. 
Therefore we obtain: 


(5.35) 


(5.36) 


d . dT . dT A .. dr, 

dx^dqj dq k 9q k 


(5.37) 


Applying Newton’s law to the right hand side and introducing the generalized forces 
of 5.28 and 5.29 


%r Q '' +Q ' 


We arrive at Lagrange’s equations: 


d_ 

dt 


dT dT ( £ ) q (r) 

(qT" ) "n Vfc “•'Vic 

dq k dq k 


k = 1,2, . . . , n 


(5.38) 


(5.39) 
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e generalized reaction forces Q[ can be eliminated if the given mechanical sys- 

t ! m T C l OnfigUratl0n Can be defined in terms of independent generalized coordinates 
q r This is only possible if the kinematical constraint equations can be expressed in 

holonomic form as in Equation 5.11. Using the principle of virtual work of Equa- 
tion 5.19 we have: 

N n — m 

m = E R. <r. = y Qf ) 6 = 0 ( 5 . 40 ) 

1=1 i=i 

where m is the number of holonomic constraint equations. Now the 6q * can be 
chosen independently and the generalized reaction force associated with each gen- 
eralized coordinate must be zero. This is so because we can now let all 6 q m be zero 
except for any one 6 q* which will be chosen not to be zero. 

For holonomic systems, the Lagrange equations assume the form: 

d ( dT \ dT ^ 

dt^dqj dq k ~ * = 1 > 2 »---.n (5.41) 

In this case, the number of generalized coordinates is equal to the number of degrees 

of freedom of the system. The superscript E has been dropped for simplicity of 
notation. 


Usually some or all external forces can be derived from a potential energv I 
such that 


F = -gradF (5 .42) 

where V = F(r x , r 2 , . . . , r„, i ) = V ( 9l , ? 2> . . . , q n , t). 

We calculate the generalized forces arising from such an irrotational force field 
lrom the associated virtual work: 


N 


6W 


N 


- 52 F * . = - 52 S rad Vi ■ 6r, 


1 = 1 


N 


i = l 
n 


= -52 grad 1] • (j; jp- 6 q k ) 

k=i d ^ 


1 = 1 


= - £ (52 ^ ad Vi ' dTi " c ^ dv 


k= 1 i'=l 


9q k 


) % = - 52 -kj- % 

k = i dqk 
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The generalized force is then: 


NOTE: 



(5.43) 


If the potential energy V is independent of time, the force field is called con- 
servative. This does not imply that the mechanical system is conservative. In this 
case, the Lagrange equations are written in the form: 


±(?L,dT d V 
dt dq k dq k dq k 

In theoretical mechanics, it is often convenient to 


k — 1, 2, . . . , n 
introduce the Lagrangian 


(5.44) 

function 


L as follow s: 


L = T-Y (5-45) 

If furthermore all forces are derivable from a potential energy, Lagrange equations 
become 

= 0 k = 1 , 2 , . . . , n (5.46) 

dt K dq k dq k 

This is the standard form of Lagrange’s equations for holonomic systems. 


NOTE: 


For practical applications of the form of Equation 5.44 is more useful. 

In some electromechanical systems, the potential function I can also be de- 
pendent on the velocities of the generalized coordinates, i.e., V = \ (qj,qj, 0 - The 
generalized forces are then obtained bj the prescription. 


dv d ar 

Qk ~ ~di k + dt [ dq k } 

It is apparent that in such a case, Equation 5.46 is still applicable. 


(5.47) 


Examples 

a) Scleronomic System 

Consider the double pendulum swinging in the vertical plane where the two 
masses mi and m 2 are connected by massless bars of length and ^ 2 as shown in 
the figure. 
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m 3 g 

Figure 5.4: Double Pendulum 

The configuration of this system can be defined by the Cartesian coordinates 
(x u yi) of mass mi and ( x 2 ,y 2 ) of mass m 2 . These four coordinates must satisfy 
the kinematical constraint conditions. 

(“) *\ + y\ = 

(X 2 - Xyf + ( 3/2 ~ Vlf = l\ 

The system has therefore two degrees of freedom. As a consequence, it can be 
defined by two independent generalized coordinates for which we choose the two 
angles 9 X and 9 2 which the two bars make with the local vertical. These two angles 
are related to the four rectangular coordinates as follows: 

W *1 = l\ sin 9 1 yi = lx cos 9 { 

x 2 = ^2 sin 9\ T -^2 sin $2 

y 2 = t 2 cos 9\ +l 2 cos 9 2 

The kinetic energy of the system is: 

( c ) T =\ m i (*? + V 2 i) +\ m 2 (*2 + yl) 

This can be expressed in terms of the generalized coordinates by differentiating Eqs. 
(b) with respect to time which yields after some algebraic manipulation: 

( d ) T= l - m^lj 1) 2 + \ m 2 [(tjj* + (£ 2 9 2 f 

+2 lx l 2 9x 9 2 cos{9 2 - tfi)] 
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The potential energy of the system is: 

(e) V = -m 1 gyi - m 2 gy 2 

Expressed in terms of generalized coordinates it is given by: 

(/) V = — migti cos0j — m 2 g{t\ cos0j + £ 2 cos0 2 ) 


NOTE: 

The potential energy was referenced to the level y = 0 at which V' = 0. 

The equations of motion are now obtained by performing the differentiations 
called for in Equation 5.44. The right hand side of Equation 5.44 is zero in this case 
because the external force was expressed by its potential energy. 

flt-equation: 

(fir) - fir 'nr ~ ( mi + m2 )^! + ^ cos (^ 2 _ 

dt o@i Ov\ Ou i 

-77i 2 £\ 62 sin (0 2 - #i) + (mi + m 2 ) l x g sin 6 1 

= 0 

02-equation: 

d , dT . dT d\ - 2 /i‘ c c 77 (a d \ 

17 (tt) - it + = m 2 + m 2 e 1 e 2 e 1 cos( 0 2 - 0 X ) 

al Ou 2 Ou 2 Ou 2 

• 2 

+m 2 £\ £ 2 01 sin (02 — ^1) -h m 2 ^2 P sin 0 2 
= 0 

b) Rheonomic System 

As an example of the rheonomic constraint, i.e., a constraint which contains 
time explicitly, consider a simple pendulum of mass m attached to a string whose 
free length £ can be varied by pulling the end A. 

The position of the mass is defined by the length £ and the angle 0. But since £ 
is a function of time depending on the motion of A, the system has only one degree 
of freedom with 0 being the only generalized coordinate. 
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Figure 5.5: Rheonomic Pendulum 
The constraint equation is 

x 2 + y 2 = £ 2 (t) (5.48) 

The kinetic energy of the system is: 

T = | m {i 2 + y 2 ) = l - m (P + i 2 9 2 ) (5.49) 

The potential energy of the system is: 

V — — mg y = —mg t(t) cos 9 (5.50) 

Substituting these terms in the Lagrange Equation 5.44 we obtain 

= m t 2 6 + 2 mi £9 + mg £ sin 9 = 0 (5.51) 

It is interesting to introduce a small angle approximation such that sin 9 « 9 and 
write the equation of motion as 

9 + j0 + 9 - 9 = ° (5.52) 

Comparing this result with the well-known equation of a damped linear system, 
it is seen that the second term containing the rate of change of the pendulum length 
acts as an effective damping. In fact, the damping will be positive for l > 0 and 
negative for l < 0. It is obvious that the mechanical energy of this system is no 
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longer conserved. The agent producing this change in energy is the centrifugal force 
in the string. 

Another example of a rheonomic system is a launch vehicle whose attitude is 
controlled by proper thrust vectoring. We consider the motion of the vehicle in an 
inertial reference plane X-Y, called the yaw plane. The vehicle is assumed to be 
rigid and only subject to the thrust force of a single gimballed engine. 



Figure 5.6: Attitude Control 

The system is defined by the position of its mass center relative to the inertial 
X-Y plane, its yaw angle \ p and the engine deflection 6 C . However, the engine 
deflection is not a generalized coordinate because it is controlled by the autopilot 
which operates from signals generated by position and rate gyros located on the 
vehicle. It is therefore a rheonomic constraint and referred to as “control variable.' 
Consequently, the system has three degrees of freedom. Its kinetic energy is given 
by: 

T = l - m (x 2 + y 2 ) + | I ip 2 (5.53) 

where I is the moment of inertia of the vehicle about its mass center. 

The generalized forces associated with the three generalized coordinates are: 

Q x -F cos (ip - 6 C ) , Q y = F sin (ip - 6 C ) , Q w = F £ sin 6 C (5.5-4) 

where £ is the distance of the mass center from the engine survival point. 
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From the Lagrange Equations 5.44 we obtain the equations of motion as (con- 
sider m to be constant): 


m i = F cos (tp — 6 C ) 
m y = F sin (^ - <5 C ) 

I ij = F £ sin 8 C 

For preliminary control dynamics analysis, it is useful to assume small angle 
deflections. A very simple and effective altitude control can be obtained by a so- 
called attitude/attitude rate control law expressed as: 


6c = -{h ip + xp) 

where 


k^ — position gyro gain 
k 0 = rate gyro gain 

Assuming small engine deflection angle 8 C the rotational (attitude) equation 
becomes by introducing the above control laws: 


/ i> + {F £ kj i> + (F l **) rp = 0 

This is seen to be the differential equation for a damped oscillation. It should, of 
course, be realized that this is the simplest mathematical model and its sole purpose 
was to expose the rheonomic nature of a feedback control system. 


5.5 Lagrange Equations With Reaction Forces 


If the system is nonholonomic, a reduction of the generalized coordinates to inde- 
pendent generalized coordinates is not possible and one has to use the Lagrange 
equations in the form given by Equation 5.39. This form of the Lagrange equations 
allows furthermore to calculate internal dynamic stresses of critical system elements. 
Sometimes it is also possible to introduce surplus coordinates in order to simplify 
the equations of motion. This procedure provides a means to trade off simplicity 
versus number of equations to be solved. 
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The Lagrange multiplier method provides an elegant and efficient way to solve for 
these unknown reaction forces caused by the kinematical constraints. The method 
is equally applicable to holonomic and nonholonomic systems. However, if applied 
to holonomic constraints as given in Equation 5.11 they have to be written in 
differential form as virtual constraints: 

6< t>j = Y, IT" 8 9k = ° j ~ 1,2, ...,m (5.55) 

k= i a 9k 

Recall that virtual displacements are taken irrespective of time such that 8t = 0. 

According to the principle of virtual work the mathematical problem is that 
of determining the stationary value of a function namely the virtual work, of the 
reaction forces subject to the kinematical constraints of Equation 5.55 


= 0 (5.56) 

k = ] 

In the subsequent analysis, we can treat holonomic and nonholonomic systems 
alike. The only difference between the two systems is that for holonomic constraints, 
the constraint coefficients c j k take the form of partial derivatives such that: 


_ d+j 
c,k ~ % 

Applying the Lagrange multiplier rule, the stationary value of the virtual work can 
be determined by multiplying the m virtual constraint equations 

n 

Y = 0 j = l,2,...,m (5.57) 

k- 1 

by the m Lagrange multipliers A j and add them to Equation 5.56 to obtain: 

5Z(Qfc fl ^ + Ai c lk + A 2 c 2 /t + . . . A m c mfc )6g* = 0 (5.58) 

*=i 

In this sum of n terms, we can now select the m Lagrange multipliers in such 
a way that the last m terms vanish. The remaining (n-m) terms contain then only 
(n-m) variations 6q k which are independent. Therefore, each associated coefficient 
has to be zero. 

Thus, we obtain the set of equations for the reaction forces: 
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(5.59) 


m 

Qk R) = ~T, c >k ^ = 1,2, ... ,n 
j=i 

NOTE; 

Some authors multiply the constraint equations ( 5.57) by the Langrange multi- 
pliers and subtract them from the virtual work equation ( 5.56). This only changes 
the sign of the right-hand side of Equation 5.59. 

Despite the fact that these n equations all have the same form, it is well to 
remember, that they have different origins. The last m equations hold because we 
selected the m Lagrange multipliers to make them true. The first (n-m) equations 
are true because the Lagrange multipliers are selected to make the associated virtual 
displacements independent. (See Appendix for details on Lagrange Multiplier Rule.) 


We have now (n+m) unknowns, namely the n generalized coordinates and the 
m Lagrange multipliers. But we have also the same number of equations, namely 
the n equations of motion from Equation 5.39 and m equations of constraint from 
Equation 5.57. 

For holonomic systems, it is useful to express the constraint equations in terms 
of velocities as: 


A ,d<t>j . d<j>; n 

£ + aT = 0 3 = 1,2 m (5 ' 60) 

In special cases the set of equations can be solved by ad hoc elimination and substi- 
tution methods. However, a systematic algorithm can be setup taking the following 
steps: 

1) The Lagrange equations are a set of n second order ordinary differential equations 
which can be written in matrix form as: 

A/x = Q ,4 + Qy -|- (5.61) 

where M is the (n x n) generalized mass matrix which is nonsingular for a proper 
mechanical system. The vector x is the (n x 1) column matrix of the generalized 
coordinates. The three terms Q,j,Q/ and represent the generalized applied 
force, the generalized inertia force and the generalized reaction force, respectively. 
All these forces depend only on x, x and time. 
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2) Introduce the constraint equations (holonomic or nonholonomic) in matrix form 
as: 


C x + b(t) = 0 (5.62) 

where C is the (m x n) constraint matrix. The generalized reaction force of Equa- 
tion 5.59 is likewise expressed in matrix rotation as: 

Qfl = -C t A (5.63) 

where A is the (m x 1) column matrix consisting of the Lagrange multipliers. 

3) Differentiate the constraint Equation 5.62 to obtain: 

C x + Cx + b(() = 0 (5.64) 


4) Premultiply the equations of motion ( 0.61) by the inverse generalized mass 
matrix M~ x and substitute in Equation 5.64. This yields: 

CM- l (Q A + Qi + QR) = -Ck-b{t ) (5.65) 


5) Using Equation 5.63 we solve the preceding equation for the Lagrange multipliers: 

A = (CM~ l C T )~ l {Cx + b(<) + CM~ X {Q A + Q/)} (5.66) 

The reaction can now be readily obtained by going back again to Equation 5.63: 


Q fi = -C T {C {CM~ 1 (Qa + Q /) + Cx + b(<)} (5.67) 


Examples 

a) Nonholonomic System 

A uniform sphere of mass m and radius a rolls without slipping on a plane 
horizontal surface. 



i 



(5.68) 


The kinetic energy of the Sphere is given by: 

T = ^m(x 2 + y 2 + z 2 ) + ^ Iu) 2 

The motion of the center C of the sphere is constrained by the roll condition 

v c = U) X R 0 where R 0 = aj (5.69) 

The rotational velocity u) of the sphere has to be expressed in components of the 
inertial reference frame x, y, z. 

We can relate the unit vectors ei,e 2 ,e 3 of a body-fixed system to the unit vectors 
i, j, k of an inertial frame and obtain the following result: 


UJ — + U>2 e 2 “h ^3^3 

(^1 All "I" ^2-^21 + w 3^3i)i + + U 2 A 22 + ^3^432)j 

+(wi-4i3 -f W2T23 + W 3 T 3 s)k 

Likewise, the velocity v c of the sphere’s center is expressed in inertial components 
as: 

v c = ii + yj + ik (5.70) 

Substituting the above two equations, Equation 5.69 yields: 


x — — a(wiTi 3 + u; 2 T 2 3 + W 3 . 433 ) 


y = 0 


(5.71) 


i — a(uJiAii + UJ 2 A 21 + ^ 3 ^ 431 ) 

The middle equation represents a holonomic constraint which states geometri- 
cally that the center of the sphere maintains a constant distance above the horizon- 
tal surface. It can be easily taken into account by just setting y = 0 in the kinetic 
energy expression of Equation 5.68. 

The other two constraint conditions can be expressed in terms of generalized 
velocities by introducing the classical or modern Euler angle system. Using tie 
modern Euler angles we obtain: 
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i 4 -a ,4 13 <j) 4- (sin 9A i3 — cos 9 sin <f>A 2 3 - cos 9 cos <f>A 33 ) 


4-(sin <j)A 33 — cos^».4 23 )0 = 0 
i 4- (— A n )<j> + (sin 6 A u - cos<^4 2 i — cos 6 cos <f>A 31 )ip 
4-(sin ^L4 31 — cos <j) 4 2 i)0 = 0 

This set of constraint equations can be expressed in matrix form as: 

C q = 0 C = (2 x 5)malrix (5-72) 

where = [x z $ ip 9}. 

The equations of motion are then obtained by expressing the kinetic energy in 
terms of the generalized coordinates 

T = ^ m (z 2 4- i 2 ) 4- ^/(0 2 4 - ip 7 + 4> 2 - 2<j>ip sin 0) (5.73) 

and performing the differentiations required by the Lagrange equations. 

m x = — (c u Ai 4- c 21 A 2 ) 
m y = — (ci 2 Ai 4- c 22 A 2 ) 

1(9 + (pipcos9 ) = — (c 13 Ai 4- c 23 A 2 ) 

I(<j> — -0 sin 9) — ip9 cos 9) = — (c 14 Ai 4- c 24 A 2 ) 

I (ip — 0 sin 9 — ip9 cosf?) = — (c 15 A : 4- c 25 A 2 ) 

The solution of this set of equations could formally proceed along the steps 
outlined above. The initial conditions would have to be chosen to be compatible 
with the requirement of pure rolling without slipping. 

b) Holonomic System (Scleronomic) 

Two particles tji\ and m 2 are connected by a massless rod. They move in a 

vertical plane under a frictionless constraint which keeps m t on the horizontal x- 

axis and m 2 on the vertical y-axis. Calculate the stress in the rod. 
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Figure 

5.7: Constrained Dumbbell 


The kinetic energy is: 



and the potential energy: 

T = i m (f6 2 + l 2 ) 

(5.74) 

V — —mgt cos 8 

The kinematical constraint condition expressed in terms of the velocity: 

(5.75) 


£ = 0 

(5.76) 


The Lagrangian equations yield the following two equations of motion: 


l 0 + g sin 0 = 0 
m l — m 16 2 — mg cos 6 + A = 0 

The first of these equations can be solved by itself and furnishes the time histories 
of the angle 6 and the angular velocity 6. With these known, one can solve for the 
reaction force associated with the constraint ( = £ = 0 : 
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R ( = —A = —m £ O 2 — mg cos 6 

The reaction force is negative, indicating that the rod is in tension. Physically 
interpreted it is observed that the reaction force is caused by the gravity component 
along the rod and the centrifugal face arising from its rotation. 

c) Holonomic System (Rheonomic) 



Figure 5.8: Sliding Mass 

A rigid wire is pivoted at one end so that it can be rotated in a horizontal plane 
with prescribed angular velocity u = f(t)- A particle of mass m can slide without 
friction along the wire. Calculate the reaction force of the wire on the mass. 

The kinetic energy is: 

T = i m (f 2 + r <^ 2 ) (5.77) 

The rheonomic constraint condition: 

uj = <}> = f(t) or 4>- }{t) = 0 (5.78) 

The Lagrange equations are: 


mr — mrcf) 2 — 0 
m r 2 <f> + 2 mr r <f> + A = 0 

The generalized reaction force associated with the generalized coordinate <j> is the 
reaction torque about the pivot point 0 and given by: 
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( R) 

Q<p — = mr 2 <f) + 2mrr<f) 

The second term is the Coriolis torque exerted by the wire on the sliding 


( 

mass. 
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Chapter 6 

Modal Synthesis Technique 


6.1 Boltzmann-Hamel Equations 


When applied to complex dynamic configurations, the classical Lagrange equations 
become formidably lengthy and their computer coding leads to low computational 
efficiency. This is due to the fact that they contain complete information concerning 
the dynamics and kinematics of the system. It is possible to separate these two 
aspects of a dynamical system and obtain a substantial reduction in complexity 
of the equations of motion by transforming the classical Lagrange equations from 
an inertial reference frame to a moving reference frame. The accompanying loss of 
information concerning position and orientation can be supplied by the appropriate 
kinematical differential equations governing the relation, between inertial and non- 
inertial (nonholonomic) velocities. To achieve the desired transformation we write 
the classical Lagrange equations in matrix form as: 

d ,5T , dT ^ _ 

dt W 0q “ Q_C A ~ Q + Qfl ( 61 ) 

where the partial derivatives represent the (n x 1) column matrices containing as 
elements the partial derivatives with respect to the n generalized coordinates a (i = 
1,2 ,...,n). 

The transformation from inertial velocity components to non-inertial velocity 
components is effected by the kinematical differential equations. 

f2 = .4(q)q (6.2) 

The nonholonomic velocity vector fi is, in general, composed of both linear and 

177 


\l\t 


wwwoMAirr *aa* 


PRECEDING PAGE BLANK NOT FILMED 


angular velocity components. The kinematical conditions for the angular velocities 
are given by Euler’s kinematical differential equations of section IV. The transforma- 
tion from translational inertial velocities to non-inertial velocities is accomplished 
by the direction cosine matrix A such that: 

v = .4x (6-3) 

where x T = ( x,y,z ) and v T = (t^ v 2 u 3 )- The kinetic energy of the system when 
expressed in nonholonomic velocity components is denoted by 


T’(S2, q) = T(q, q) (6.4) 

The superscript star should indicate that the mathematical form of the transformed 
kinetic energy differs from the original one although its scalar value is, of course, 
the same. 


The differentiations required by the Lagrange equations are first given in scalar 
form as: 


^L = y = y — Aik (6.5) 

dq k y dQ ( dq k V Ml 

where A ik are the components of the kinematical transformation matrix -4(q) of 
Equation 6.2. Also: 


ar_^ or an< dr 

dq k f dQi dq k + dq k 
Equation 6.5 can be directly written in matrix form as: 

dr _ t dT ‘ 
a^ _ ■ on 

To be able to express Equation 6.6 in matrix form we introduce the Jacobian matrix 


( 6 . 6 ) 


(6.7) 


J = 


' gQi 

an t 

an, * 

dq\ 

dq? 

dqi 

dSlj 

dUj 

dUj 

_ dq ) 

dq? 

dqi . 


dfl 

dq 


Equation 6.6 can then be written in matrix form: 

dT T dT’ dT' 

aq df2 + dq 

Substituting Equation 6.7 and 6.9 into Equation 6.1 yields: 


( 6 . 8 ) 
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(6.9) 


a t !("!)+ 4 

dv on )+ dn J dn 


dT m 

dq 


= Q-C r A = Q + Q* 


Next we premultiply this equation by (.4 _1 ) r and obtain: 


|(§£) + [(.4-,).4-F 


8T -(A~‘) T ~=(A~‘) T Q -(CA-'fX 


dn 


dq 


This is the desired transformation of the classical Lagrange equations from 
ertial frame to a moving frame. These equations are known as Bolt/mann- 
equations or Lagrange equations for quasi-coordinates. 


( 6 . 10 ) 


an in- 
Hamel 


It remains to transform the generalized forces from the inertial reference frame to 
the moving reference frame. This is accomplished by the so called quasi-coordinates. 
These serve only as a means to an end and play only a brief cameo role after which 
they disappear from the scene. 


The quasi-coordinates are defined such that their time derivatives are equal to 
the nonholonomic velocities, i.e., 



i = quasi-coordinate vector 


( 6 . 11 ) 


The virtual displacements of the generalized coordinates q { can then be related to 
the virtual displacements of the quasi-coordinates & by Equation 6.2 as: 


6£ = A(q)6q (6.12) 

Physically interpreted the virtual displacements of the quasi-coordinates are in- 
finitesimal translations along or rotations about the instantaneous axes of the mov- 
ing reference frame. Calculating now the virtual work of the external forces Q 
through a virtual displacement we obtain: 


SW = (6q) r Q = (A~'6t) T Q = (^) r (^- 1 ) r Q = (^ T )k (6.13) 

where k is the generalized force associated with the quasi-coordinate £. Mathemat- 
ically expressed: 


k = (A- l fQ = (A T rq (6.14) 

The generalized reaction force in the moving reference frame can be obtained by 
transformation of the constraint condition. 
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(6.15) 


Cq + b(<) = 0 

Using the kinematical differential equation ( 6.2) we obtain: 

C(A~ l n) + b{t) = Bfl + b(i) =0 (6.16) 

where B = CA~ l is simply the constraint matrix relating the nonholonomic veloci- 
ties of the system elements as required by the kinematical constraints. 

The Boltzmann-Hamel equations are therefore finally: 


— ( 211 ) + {(A - J)A~ l ] T — 

dt K df 2 ^ 


- U ' 1 ) 7 




dr 

dq 


k — = k H- k ji 


(6.17) 


Judged by their outward appearance the transformed equations ( 6.17) seem to 
be much more complicated than their classical counterpart given in Equation 6.1. 
However, their hidden simplicity will be brought to light when the detailed steps of 
introducing the above-mentioned moving reference frame are carried out. To this 
end we introduce the linear and angular velocity of the moving reference frame as: 



’ Vi 


" Wi 

V = 

v 2 

= 

U) 2 


. ^3 - 


UJs 


The column matrix Q can be written in partitioned form as: 


n = 


U) 

V 


It can be shown by some rather lengthy but straight-forward process that the fol- 
lowing relationship holds: 


where 


lj = 


[(.4-J).4-f = 


0 

^3 

U>2 

U^3 

0 


— U>2 

Wl 

0 


UJ j 

V 

0 1 

u 


■ 0 - 

V = 

^3 


-v 7 


0 

Vi 


v 2 

-Vi 

0 
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Next we assume without any practical loss of generality that the kinetic energy of 
the system can be expressed as: 


T = T(f2,ri,ri) (6.18) 

where 7} are generalized coordinates defining the system configuration relative to 
the moving frame. The superscript star has been dropped for simplicity. With this 
assumption the term dT* jd q in Equation 6.17 is equal to zero. The equations of 
motion can now be written as two vector equations: 

A) Translational Equations: 

d.dT. dT „ „ x 

5? ( 57 ) + u ’ x a7 = F + F '' (619) 

B) Rotational Equations: 


d .dT. dT 

* ( ^ ) + “ x ^ + v 


dT 

x — — L + Lfl 

dv 


(6.20) 


The right-hand side of Equation 6.19 are the forces acting on the system in the 
direction of the instantaneous moving frame axes, that of Equation 6.20 are the 
moments of these forces about these instantaneous axes. The latter can be derived 
easily by calculating the virtual work done by the forces through the virtual rota- 
tions of the moving frame. Using the quasi-coordinate definition of Equation 6.11 
the virtual displacement of a point located at R in the moving frame is given by: 


8xft = 6£ r x R Remember: v = u; x R 
The virtual work is therefore: 


( 6 . 21 ) 


8W = F • 6x r = F • (d£ R x R) = (R x F) • d£ R = L • d£ R (6.22) 

mm 

The moving reference frame does not have to be a body-fixed frame but can be 
any conveniently chosen frame relative to which the system motion is to be defined 
(“Floating” frame). 

The Lagrange equations for the relative motion retain their classical form. Since 
the relative motion is often associated with flexible system components we call them 
flexibility equations: 
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C) Flexibility Equations: 


d_ dT _ dT 
dt di] dr) 


Q + Q* 


(6.23) 


The real physical meaning of the equations of motion given by Equation 6.19, 
6 . 20 , and 6.23 can be revealed by introducing the explicit expression of the kinetic 
energy of the system in terms of the linear velocity v of the origin of the reference 
frame, its angular velocity w against inertial space and the relative motion R of the 
system masses as viewed from the moving frame. With these terms we obtain: 


T = - f [v + (u> x R) + R] dm (6-24) 

2 J 

In the subsequent mathematical manipulations the following relations and identities 
are used: 


(a) If y = (a x b) 2 then fj = 2b x (a x b). 

(/?) If y = (a • b) then f* = b 

( 7 ) u x (R x v) + v x (u> x R) = R x (w x v) 

(d) w x [R x (w x R)] = R x [u> x (u? x R)] 

(e) w x (R x R) + R x (w x R) = R x (w x R) 

It is also important to realize that all time derivatives are taken relative to the 
moving reference frame. 


We obtain: 


dT 

du> 

J R x [v + (w x R) + R] dm 

(6.25) 

dT 

W X du ~ J 

f u x |R x[v + (wxR) + R]J dm 

(6.26) 

dT 

XX lh 

= J v x [v + (u> x R) + R]dm 

(6.27) 

-) = / 
du> } J 

R x [v + (u> x R) -f R]dm 

(6.28) 

+ 

/ R x [v + (w x R) + (u> x R) + Rjdm 
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(6.29) 


d ,8T . r , . 

^■(^ 7 ) = y ( v + (w X R) -f- (u> x R) + R]dm 

dT r 

U,X 57 = y a,X f V + ( u ’ xR ) + R ^ dm (6.30) 

Collecting all terms associated with the translational equation we obtain: 

|[v + (wxv) + wx(wxR)+wxR + 2( W xR) + R]dm = F + F* (6.31) 
For the rotational equation: 

/ R X [v + (W X v) + U X (u x R) + W X R + 2(u> x R) + R]rfm = L + L fl (6.32) 

It is observed that the terms within the square bracket taken together determine 
the absolute acceleration of a mass point relative to inertial space. The first two 
terms represent the acceleration of the origin of the moving frame expressed in 
components of the moving frame. The bracketed term is thus seen to be identical 
with Equation 1.5 of section I. The translational and rotational Equations 6.19 
and 6.20 could have been obtained also by a Newton-Euler approach. The main 
difference occurs in the physical interpretation of the right-hand side. Here they 
are seen to be generalized forces associated with the virtual displacements of the 
instantaneous reference frame. Furthermore the elimination of constraint forces can 
be very efficient by accomplished by the Lagrange multiplier method. Constraint 
forces will arise when several moving coordinate frames are introduced for various 
sj stem elements and then conjoined for dynamic simulation. 


6.2 Component Modes 


In many situations the motion of flexible components can be described by the super- 
position of appropriately chosen mode (shape) functions. The success of this method 
often referred to as component mode synthesis depends largely on the proper choice 
of these assumed mode functions. These are often selected from the natural modes 
(eigen functions) of the isolated structural component (substructure) using bound- 
ary conditions which are geometrically and dynamically resembling the actual ones. 
Other mode functions can be obtained by the static deflections of the substruc- 
ture due to unit displacements or unit forces imposed upon suitable coordinates. 


183 


Many other types of component modes have been advocated to describe the flexural 
motion of the substructure (Ref: Roy R. Craig Jr. “Structural Dynamics” 1981). 

In all these cases the position of a mass element in the moving reference frame 
can be defined as: 


R = r 4- V\( r ) */.'(<) (6.33) 

x 

where V\(r) is the vector mode function and T] ,{t) its associated generalized coordi- 
nate. The vector mode function is often specified in terms of a translational mode 
of the center of mass of a mass element and a rotational mode about its mass center 

such that 


V>,(r) = «/>,(: r) + <f>[ * p ( 6 - 34 ) 

where p is the position of a mass particle relative to the center of mass of the 
mass element. The position vector r denotes the location of the undeformed mass 
element. 


The flexibility Equations 6.23 can then be further manipulated using the follow- 
ing partial differentiations: 


dT_ 

dr]. 


dr_ <m 

dR ' dr], 


- = - / V\ • {“> x [ v + ( 

i J 

dT dT dR [ . r ,, 

= — • IT = / • [v + ( 

dr], (9R dr], J 


X 


x R) + R]| dm 
R) + R]dm 


— (^— ) = / V*, 1 [v + (<*> x R) + (w x R) + R]dm 
dtdr], J 


Collecting all terms we obtain Equation 6.23 in the form 


(6.35) 

(6.36) 

(6.37) 


J ip - ■ [v + (uj xv) + w x ( u> x R) + (u> xr) + 2(w x R) + R ]dm — Q,' + Q, (6.38) 

The bracketed term is again seen to be the total acceleration of the mass particle 
relative to inertial space. Therefore the generalized inertial forces are given by the 
scalar product of the vector mode function with the total acceleration and summing 
over all mass elements. The generalized forces are obtained via the virtual work 
through a virtual displacement of the generalized flexural coordinate rj,. Therefore: 

Qi = F • V\( r ) an d Q\ R) = Ffl • VhM (6.39) 
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Another contribution to the right-hand side arises from the strain energy of the 
flexible component. For a beam, for instance, undergoing a transverse deflection 
<P,(rmO the strain energy is given by 

V - \ J Vifdx (6.40) 

Its generalized force is obtained by partial differentiation 

fj? = Ql B> = j Eltfdx ,, (6.41) 

Similar expressions can be obtained for the strain energy of other system elements 
which undergo a deflection as dictated by the assumed mode function. 

It is often advantageous to select mode functions which are orthogonal with 
respect to the mass distribution of the component. This eliminates or reduces the 
number of dynamic and static coupling terms. 

As an example, consider the following integral which is encountered when Equa- 
tion 6.38 is expanded into its various components: 

1 = J (<f>i + <t>\ X P) ■ ; + <t>'j X p)dm (6.42) 

j 

Orthogonality of the mode functions implies that 

J ( 0 , + <t>i x p) • {<f>j + </>' x p)dm = 6 {j ( 6 . 43 ) 

where is the Kronecker symbol. 

Therefore, the integral of Equation 6.42 reduces to 

1 = J (<t>i + <t>[ x pfdm ( 6 . 44 ) 

= <t>[)dm = 1 (6.45) 

where J = p 2 £ — (pp). 

The integral I is the generalized mass of the mode function normalized to unity. 


If the systems contains continuous and discrete elements simultaneously, the 
above mathematical formulations remain unaltered when all integrals are defined 
as Stieltjes integrals. 
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6.3 Applications to Aerospace Systems 


Example 1: 

Spherical Pendulum 



The direction of the pendulum 
is defined by the two Euler angles 
0 and V’ 

/. = ml 7 + l = h 

h = I 


The rotating reference frame is defined such that the e 3 axis is along the tether. 
The rotational motion is governed by the rotational Equation 6.20 vt hich becomes 


T-ui + u}xI-u} = 1j + L fl — L B X 

The position of the pendulum in inertial coordinates is given by: 


x = l sin 6 cos xp 
■y — £ sin 0 sin V 
z = £cos^ 


Its angular velocity is given by: 
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u>i = — xjj sin 9 

u> 2 = 9 

0 J 3 = ificos9 

Because <j> = 0 then exists a kinematical constraint between the angular velocities. 
It is: 


cos 9 + u > 3 sin 9 = 0 
which leads to the constraint matrix: 

B = [cos 9 0 sin 9] 

The torque acting on the pendulum comes from the gravity force and is: 

L = r x F = £e 3 x mgk = — mg £ sin 0e 2 
The equations of motion can now be obtained in scalar form: 


IiLJi + 1 ^ 2 ^ 3(73 — I 2 ) — A cos 9 
/ 2 W 2 + — J$) = — mg£ sin 9 

73 W 3 + 1 ^ 1 ^ 2(72 — 7J = A sin 9 

To eliminate the Lagrange multiplier we multiply the first equation by sin 9 and 
the third equation by - cos# and add. We obtain the following two differential 
equations: 


(7 + m7 2 sin 2 6)ij) + 2 m£ 2 Tp 9 sin 9 cos 9 = 0 
I 2 0 — m£ 2 xjj 2 sin 9 cos 9 + mgt sin 9 = 0 

Those equations could have been obtained by using the classical Lagrange equations. 
It can be seen that for 7 << ml 2 these equations become almost singular for small 
9 and present serious computational problems. 

To avoid this singularity it is necessary to express the direction of the pendulum 
in terms of the modern Euler angles 9 and <f> (ip = 0). 
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The angular velocity of the pendulum is then 


= (j) x = £ s\n 8 cos <f) 

u >2 = 8 cos (f> and y — — £ sin <f> 

w 3 = —8 sin <p 2 = £ cos 8 cos <f> 

The constraint imposed upon the angular velocity because of = 0 is: 

w 2 sin <f> + u >3 cos <j) = 0 
which leads to the constraint equation: 

B = [0 sin <{> cos <f>] 

The torque acting on the pendulum is given by 

L = 7e 3 x mgk = -mg£ sin 0e 2 - mg£ sin <j> cos 

The rotational motions are therefore: 

I 1 uj 1 + u) 7 u) z {h ~ h) = —mgt s\n<f>cos8 
u^ 3 (-f i -^ 3 ) ~ m.g£ sin 8 T A sin 

/ 3 W 3 + Wi(x>2(f2 — ^l) = A COS (f) 

To eliminate the Lagrange multiplier we multiply the second equation by cos <j> and 
the third equation by sin <f> and obtain again two differential equations: 

(m£ 7 + I)(j> + mi 7 8 7 sin (f> cos </> = -mgl sin <f> cos 8 

(mi 7 cos <f> + 1)8 - 2 m£ 7 8<j> sin <f> cos (j> = -mgl cos <f> sin 8 

It is seen that these equations are well-behaved for small angles 8 and <j> and in fact, 
reduce to the simple pendulum equations: 

(7 + mt 7 )(j> + mg£<t> — 0 
(7 + m£ 7 )8 + mgl8 = 0 

Example 2: 
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A dumbbell is elastically coupled to a flat disk at its mass center and can rotate 
in a plane through an angle a as shown in the figure. 


x 3 



x i , x 2 , 1 3 body-fixed 
reference frame. 
m i = m? = m 

The origin is at the center of mass 
Jo = 2ml 7 


Moment of Inertia of Disk: 


J = + J9e 2 e 2 + Ce 3 e 3 

Moment of Inertia of dumbbell: 


T 0 = Jo cos 2 + /oe 2 e 2 + / 0 sin 2 ae 3 e 3 — 7 0 sin a cos aeie 3 

Location of dumbbell masses and m 2 : 
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Angular Velocity of Disk: 


Ri = t cos ae 3 + £ sin ae! 
R 2 — —{£ cos ae 3 + l sin aei ) 


Kinetic and Potential Energy of Dumbbell: 


1 


1 


T 0 = -/ 0 (d + D 2 ) + sinQ; - cosq ) 2 

To = ~ka 2 k = spring constant 


Rotational Equation: 


X • /? “I - /? x X * /? Xq * *f? "4" if? x X"o * 7? •+■ 2Rj x (/? x V[ )m ^ 

+2R 2 x (/? x v 2 )m 2 + (R[ x ai)m[ + (R 2 x a 2 )m 2 = 0 

where v and a are the time derivatives of the position vector R of the dumbbell 
masses relative to the body-fixed frame. 

Flexibility Equation: 


dV da } 

dl o 

da 

d_ on 

di da 


da + da 

— ?o(d + D 2 ) 

) = ?o(<* + f?2) 


dTo 

da 


7o(Q 3 sin a — f?! cos a)(Q 3 cos a + fii sin a) 


d [Vi 

da 


= k a 
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Information about the attitude of the main body can be supplied by the Euler 
kinematical equations. For setting up perturbation (linearized) equations about a 
nominal spin condition it is, of course, necessary to use the modern Euler angle 
system to avoid the gimbal look singularity. 

Example 3: (Slosh Model) 

A satellite contains an internal linear oscillator system located at r 0 . having a 
spring constant k and a damping constant c. The motion of the oscillator is given 
by its displacement vector r,. 


3 



m, = mass of oscillator 
M = total mass of satellite including m, 

r OJ T r. 


The body-fixed reference frame origin coincides with the center of mass of the 
satellite for r, = 0. The location of the center of mass of the satellite for a displaced 
oscillator mass is: 


Translation: 


M 


l c = J R dm = m,r. 


The translational Equation 6.31 contains the following terms: 
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v 0 + (u) x v 0 )]dm = A/(v 0 + w x v 0 ) 


J u) x (w x R )dm = u» x (u> x j Rdm) = M u? x (w x £ c ) 
^(ti; x TV) dm = oj x J Rdm = A/(u> x £ c ) 


2 



2m , (oj x f,) 


J i,dm = m,r. 


Combining these terms yields: 


M(v o + w x v 0 ) + A/w x (u) x * c ) + A/(u; x £ c ) 
+2 m, (w x r,) + m, t, = F 


Rotation: 

The rotational Equation 6.32 contains the following terms: 


J R x (v 0 + u x v 0 )dm = M[i c x (v 0 + u> x v 0 )] 


J R x [u> x (u> x R)]dm = w x T' w 


/ 


Rx (u x R)dm — T ■ uj 


2 


/ 


R x (u> x T,)dm 


= 21, x (u> x r,)m. 
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/Ex T,dm = (f, x T,)m, 

Combining these terms yields: 

I'W + «xL« + M[£ c x (v 0 + u> x v 0 )] 

+2m,£, x (u> x r,) + m,(£ , x r,) = R x F 

The effect of the acceleration of the origin of the body fixed reference frame can 
be eliminated from the above equation by substituting the equation before it into 
it: 

lu) + w x 1 • u> + £ c x [F - Mu) x (u) x £ c ) - A/(u> x £ c ) - 2m, 

(u: x r.) - m,r,] + 2 m,£, x (u> x f,) + m,(£, x r f ) = R x F 

If we define the dyadic: 

I 2 = M[£ 7 c T - (£ c £ c )} 

the moment of inertia about the center of mass of the satellite can be written as: 

Ic = I - 

The rotational motion equation is then 

I c '« + wxI c ’W + 2 m,(£, - £ c ) x (u> x r,) 

+m,{£ , - £ c ) x r, = (R - £ c ) x F 

Flexibility: 

The motion of the linear oscillator can be written in the form as: 

£. = r 0 , + e, r, e, = unit vector 

It is seen that the unit vector e, can be looked upon as the vector mode func- 
tion of the linear oscillator and its displacement r, as the associated generalized 
coordinate. 
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The flexibility equation (6.39) contains, therefore, the following terms: 


e, • [v 0 + (w x v 0 )]m. 


e,'[wx(«x t.)]m. 


e, • {u> x l,)m. 


2e, • (u> x T,)m, 


e, • r,m, 

Collecting these terms yields: 

e, • {[v 0 + (a? x v 0 )] +wx(wxf,) + (wx l,) 

+2(u> x r,) + r 4 } m, — — k r, — c r. 

Substituting the translational equation eliminates the effect of the acceleration 
of the origin: 


f r F 71/ 2m, . 

e, ■ j[— u> x(wx l c ) - (u> x l c ) - — (a> x r.) 

— ~r,] + w x (w x + (w x £,) -f 2(u > x r,) + r, j m, = — k r, — cr. 

Observing furthermore that e, ■ (w x r.) = 0 we can rewrite this equation in the 
final form: 


• {(1 - + w x (/, - l c ) + u x [w x (*. - / c )] | 


m, k c . 

= - T7 e * ‘ F r, 

M rri B m a 
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The above equations of motion representing three differential equations of first 
order and one second-order differential equation together with the appropriate kine- 
matical equations determine the dynamics of the system. 

Example 4: Rolling Coin 

A classical example of a nonholonomic system is the coin rolling on a rough 
horizontal plane. 




We choose a coordinate system 1, 2, 3 with the origin at the mass center of the 
coin. The 3-axis is the axis of symmetry; the 1-axis is in the plane of the coin and 
remains horizontal. It is important to notice that this coordinate system is not a 
body-fixed system but is “floating” relative to the body. Therefore, the angular 
velocity of the coin is different from that of the floating coordinate system. 

Angular Velocity of Coin: 


fl — + J? 0 

where f2 0 = xjje 3 (6.46) 

and u> the angular velocity of the coordinate system. 
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RqII Condition; 


v = ft x R 0 where R 0 = Re? (6.47) 

R 0 is the position vector drawn from contact point O to the mass center of the coin, 
and v is the velocity of the mass center. 

The velocity components for the roll condition are therefore: 


V\ = —f2 3 R, v 2 = 0, u 3 = S2 l R (6.48) 

For calculating the constraint forces we need the instantaneous constraints of the 
system. They are obtained by performing virtual displacements of the coordinate 
system in conformity with the kinematical constraints which apply at that instant 
of time. 


Constraint Condition: 
or in component form: 


v = u> x R 0 


(6.49) 


Vi = —W 3 R, v? = 0, u 3 = u>iR (6.50) 

Notice that the constraint condition ( 6.50) is different from the roll condition ( 6.48) 
for the coin. 


In matrix form: 


Constraint Forces: 


1 0 0 0 OR 

0 1 0 0 0 0 

0 0 1 -R 0 0 


V\ 

v 2 

u 3 

UJ3 _ 


F r = B t X 


1 

0 

0 

0 

0 

R 


0 

1 

0 

0 

0 

0 


0 

0 

1 

-R 

0 

0 




— 1 
1 


A, 1 


a 2 


i 

a 2 

. A 3 

— 

A 3 

-R A 3 
0 

J 


rtA 3 . 


(6.51) 
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Kinetic Energy: 


T=^J(\ + f2x R ) 2 dm where J R dm = 0 


Equations of Motion: 


d.dT. dT dT 

— ( — ) + wx — + v X — = L + Lfi 
dVdu> doj d\ 


Because of / TLdm = 0 


d ,dT. dT _ t _ 

dt ( fr ) + “ x a7 = r + FB 

v* f = ° 


The generalized forces L and F in the above equations are obtained by the virtual 
work done by the external gravity force when virtual translations and rotations 
about the mass center of the coin are performed. 

We find that L = 0 and F r = [0, — mg sin 8, —mg cos 8]. 

The equations of motion in component form are then: 

Translation: 


m Vi + m(cj2f3 — ^ 3 ^ 2 ) = 0 + 


m u ’2 + m{u)zV\ — ^ 1 ^ 3 ) = ~ m 9 s * n 8 + X? 


m U3 + m( CJ1V2 — u^i) — —mg cos 9 + A 3 
Rotation: (/1 = / 2 = A, / 3 = C) 


/l Qi + CiJ2^3 — — RX3 
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A fl 2 + Auj^VIi — 5^3 — 0 


( 6 . 52 ) 


C f?3 + ,4u>ifi 2 — Au>2$li = — R^3 
Combining Equations 6.51 and 6.52 we obtain: 

(/I 4- mR 2 )Qi + (C + m.R 2 )u> 2^3 — Aui^O.? — — mgRcosO 

AQ 2 + — Cwifls = 0 (6.53) 

(c + m/? 2 ) f ?3 + AujyVl? — (.4 + m/? 2 )u 2 2 fli — 0 

This equation is the moment equation about the contact point 0. 

To complete the formulation of the problem we write the angular velocity u of 
the floating coordinate system in terms of Euler angles as 

Wi = 0, ijj 7 = <f)sm6 , o> 3 = <j>cos6 

Likewise the angular velocity ft of the body as: 

Qi—8, ft 2 = <^sin0, Q 3 = <f>cos9 + ip 
A closed-form solution is only possible for certain special cases. 


C-3 
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Example 5: 


Consider two rigid bodies connected by a frictionless hinge in free space. 



We select 2 principal axes systems each located at the mass centers of each body. 
The angular velocities of each reference frame against inertial space are and ft 2 
respectively, their linear velocities are v x and v 2 . The distance from C\ to hinge 


point 0 is li and that from 0 to C 2 is l 2 . 

The equations of translation and rotation for each body are then 

Mi (vi + ft l xv l ) = F l + f[ r) (6.54) 

li • ft i T ft\ x 1\ • ft\ = lq 4- L/j ^ (6.5o) 

A/ 2 (v 2 + ft 2 x v 2 ) = F 2 + F[ ^ (6.56) 

X 2 • St 2 + ft 2 x T 2 • fl 2 — L 2 + Li 2 ^ (6.57) 

The constraint equation imposed upon the two bodies by the hinge is: 

v 2 = V! + fli x l x + fl 2 x l 2 (6.58) 


The constraint equation has to be written in matrix form to find the constraint 
matrix B. It is: 


v 2 — Tv! — All ft i — l 2 — ft 2 = 0 
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where .4 is the direction cosine matrix of body 2 relative to body 1 and t the familiar 
skew-symmetric matrix. Introducing the vector: 

fl r = [v 1 |/Ii|v 2 |/l 2 ] 

we can write the constraint equation in the desired matrix form as: 

Bn = [~A\ - AlilEifyn = 0 

The equations of motion are in matrix form: 

M n = Q (/) + Q (4) + Q (/i) 

where M is the diagonal matrix: 

M x 0 ' 

/ 1 

M 2 

0 I 2 . 

and Q^,Q^,Q^ are the associated generalized inertial, applied and reaction 

forces. The latter can be calculated by the same method outlined before. We 

obtain: 

QW = -B T (BM~ 1 B T )~ l {bm~ 1 { Q (/) + Q (y4) ) + Bn] 

The time derivative of the constraint matrix B can be obtained by noting that: 

A — —CjA 

where uj = i? 2 — -4n].4 T . It is therefore given by: 

B = [w.4|ul4^|0|0] 

N.B. 

It is important to notice that the translational Equation 6.56 is only needed to 
calculate the constraint forces and does not have to be integrated, because the 
translational velocity v 2 can be directly obtained from the constraint Equation 6.58. 
Therefore it is only necessary to integrate as many equations of motion as there are 
degrees of freedom, namely nine. 
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Example 6: 


A spinning satellite considered to be a rigid body has flexible antennas attached 
to it. 


e 2 


► e x 



The antennas are treated as uniform cantilever (clamped-free) beams. The dif- 
ferential equation of a uniform beam is: 


d i w d 2 w 

El — — r + m — — = 0 


dx i 


dt 2 


where El is the flexural stiffness and m the mass per unit length of the beam. At 
the clamped end x = 0 we have the geometric boundary conditions: 


tu(0) = w'(0) = 0 

At the free end £, we have the dynamical (natural) boundary conditions of 
vanishing moment and shear force. This means that: 


w"(l) = w"\l) = 0 

The differential equation and its associated boundary conditions furnish the 
normal modes 4> m {x) to be used in describing the total deflection of the antennas as a 
superposition of these mode shapes each multiplied by a generalized coordinate. The 
transverse deflections of each antenna are expressed relative to a body fixed frame 
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(d, e 2 , e 3 ) whose origin coincides with the attachment point A of the antenna. 
The position of a mass element of the antenna is defined as: 


t f = + ye 2 + ze 3 

According to the model synthesis method, the “in-plane” deflection y and the “out- 
of-plane” deflection z are given in the series form: 


V = Y Yi(x)Vi (0 
2 = £Z fc (ar)£ fc (0 

To define the location of the antenna mass element relative to the rigid body 
satellite we select a satellite-fixed reference frame (i,j,k) with origin at the mass 
center of the undeformed body. For the subsequent discussion it is assumed that 
the shift of mass center location caused by the antenna deflections is neglible. As 
a consequence the translational and rotational motion become dynamically uncou- 
pled. To evaluate the integrals of the equations of motion we need the following 
quantities: 


R f = l A + t f location of antenna mass element 
v f = ye 2 + ie 3 = Y Yi(x)rj t {t)e 2 + Y £*(*)& (0 e a 
a F = ye 2 + ze 3 = Y Yi(x)iji{t)e 2 + Y Z fc (ar)^(f)e 3 
= R 0 + (*T2 * £a) + ^2 x (/? x 

where = acceleration of attachment point A. R 0 = acceleration of satellite frame 
origin. 

The following equations of motion can then be established: 

Rotation: 

|Rx[17xR + flx(flxR) + 2 (f2 x v F ) + a F ]dm = L 
or 
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1- n+n xin + 2 J x (/? x \ F )dm F + J (R f x a F )dm F = L 

where dm F is a mass element of the flexible antenna. 

Flexibility: 


■ J h(x)[R A + (/? x R/r) + /? x (/? x R^) + 2(f2 x Vf) + a F ]dm F = Q t 

e 3 • j Z h (x)[R A + (J? x R f ) + n x (n x R F ) + 2 (12 x v f ) + a f ]rfm F = Q k 
The evaluation of the last two integrals is greatly simplified by observing that: 


>;•(*) = z { (x) 


J Yi(x)Yj(x)dm F = 6(j orthogonality 


j Z,(x)Zj(x)dm F = 6ij 

The acceleration Rg of the satellite origin is usually caused by thruster firings 
performed for attitude control or spin-up/spin down maneuvers. 
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Appendix A 

The Lagrange Multiplier Method 


The Lagrange Multiplier Rule finds application in determining extremes of a func- 
tion of several variables subject to constraints. 

a) Special Case: 

f = f(x,y,z) (A.l) 

Constraint: <f>(x,y,z ) = 0 
Necessary Condition for Extremum: 

df = f x dx + f v dy + f t dz = 0 (A. 2) 

From the constraint equation, we obtain: 

d<j> = <f> x dx + 4> y dy + <f> z dz = 0 (A. 3) 

Elimination of dz yields: 

dz = ~{^r)dx - {j-)dy (A.4) 

Yz *Pz 

Inserting Equation A.4 into Equation A. 2 leads to: 
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df = f x dx + f y dy - j-(<j> x dx + <j> y dy) = ( f x - )dx + ( f y - ~ )dy = 0 (A.5) 

<Pz <Pz <Pz 

Since dx and dy are independent: 

f. -&)*.=<> (A.6) 

Yz *Pz 

Define the Lagrange Multiplier A as: 

A = -f x /<f> 2 # 0) (A. 7) 

Then we obtain as necessary conditions for extrema: 


fx + Mx — 0 / t + A <j> y — 0 f z + \<p z — 0 (A. 8) 

Despite their identical outward appearance, it is important to realize that their 
origin is quite different. The last equation holds because we have selected A to 
make it true, whereas the first two equations hold because of the independence of 
the associated variables x and y. 

The Lagrange Multiplier Rule arrives at the same result by introducing an 
augmented function /' such that 


/*(*, y, z) = /(*, y, z) + \<j>{x, y , z) 

The necessary conditions for the extrema of a function subject to constraints can 
then be formulated as: 


/; = o, /: = o, /; = o 


b) General Case: 

The foregoing technique can be readily extended to the general case of n vari- 
ables. 
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m[rt 


f = f(x 

= <f>m( * 1 , ^ 2 > ■ ■ • , X n ) 

Define the augmented function 


Necessary Conditions: 


/* — / + 5Z -W« 

«'=i 


fl/‘ 

9a:* 


a** £ 


Using Matrix Notation: 


*= 1 , 2 ,..., 


Define: A r = [A, , A 2 , . . . , A m ] 
Jacobian; 


J = 


d(f> i/dxi d(j) i/dx?... d<f> \/dx n 

d<!> 2/ dx x d<j> 2 / dx 2 . . . d<t> 2 /dx n 


L dcfrrn/dx ! d(j) m /dx2 . . . d<j>m/dXm 


or 

dx 


df T 
~dx ^ ^ ^ ~ 0 


There are n + m unknowns: 


Rank J — rn 


®1 7 *^2 5 * ■ • 5 

An A 2 , . . . , A m 


There are n + m equations: 
n necessary conditions 
m constraint conditions 
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The advantage of the Lagrange Multiplier Rule is the symmetry of the formula- 
tion and avoiding the awkward elimination process of the variables. 

Example 1: 

Find the dimensions of a cylindrical can of maximum volume for a given surface 
area : 

a) Elimination Method: 

V = 7 T r 2 h S = 2ir{r 7 + rh ) = CONSTANT 


d\ = 7r(2 rhdr + r 7 dh) = 0 


(A.9) 


dS = 7ir(2rdr + hdr + rdh) = 0 = 2*r[(2r + h)dr + rdh] = 0 (A. 10) 

Eliminate dir. dh — — ^ dr 

Inserting in Equation A. 10: (2r + h)dr — 2 hdr = 0 h = 2 r 
b) Lagrange Multiplier Rule: 


V ’ = r 2 h + A (r 2 + rh) 


dV 


dr 


= 2 rh + A(2r + h) = 0 


(All) 


dv‘ 

dh 


= r 2 + Ar = 0— >A = — r 


(A. 12) 


Inserting A in Equation A. 11: 


2 rh - r(2r + /i) = 0 -» rh = 2r 2 A = 2r 
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Example 2: 


Find the dimensions of a rectangular box without a top, of maximum capacity, 
whose surface is s = 108 m 2 . 


f{x, y, z) = xyz 

<j)(x, y, z) = xy + 2 xz + 2 yz — 108 = 0 
/* = xyz + A (xy + 2 xz 4- 2 yz — 108) 


— y z + Hy + 2^) — 0 

(A.13) 

d f“ 

— — = xz 4- A(a: + 2 z) - 0 
dy 

(A.14) 

~Ir~ ~ x y + A(2x + 2y) = 0 
Oz 

(A.15) 


Multiply Equation A. 13 by x, Equation A. 14 by y and Equation A. 15 by z and add: 


x (xy + 2 xz + 2 yz) + -xyz = 0 


Using the constraint equation we find: 


. 3 , xyz 

108A + -xyz = 0 — * A = 

Inserting the Lagrange Multiplier A in ( A. 13), ( A. 14), and ( A. 15) 


1 “ 72^ + 2z ) = 0 ■"* x = y 

1 - yt(® + 2 z) = 0 


1 - —(2a: + 2y) = 0 -» 2 = 


18 


209 


x = 6 y = 6 2 = 3 

Example 3: 

Find the minimum distance of a point in a plane from the origin. 


F(x,y, z) = x 7 + y 7 + z 7 Distance 

(A.16) 

C(x, V > z ) = Ax + By + Cz + D = 0 Plane 

(A. 17) 

Augmented Function: F‘ = F + 2X0 


OF’ n 


— = 2x + 2XA = 0 

ox 

(A. 18) 

dF- 


« -2y + 2XB-0 

0 y 

(A. 19) 

dF' 


a - 2z + 2XC = 0 

Oz 

(A. 20) 

Multiply Equations A. 18. A. 19 and A. 20 by A, B ; C respectively. 

and add: 

X = — D/(A 7 + B 7 + C 7 ) 

(A. 21) 

Insert in Equations A. 18, A. 19 and A. 20: 


x M = - AD /(A 7 + B 7 + C 2 ) 

(A. 22) 

Vm = -BD/(A 2 + B 7 + C 7 ) 

(A. 23) 

z&t = -CD /{A 7 + B 7 + C 7 ) 

(A. 24) 

For minimum distance insert Equations A. 22, A. 23 and A. 2-1 in Equation A. 16: 
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Example 4: 


d M = D/y/{A* + £ 2 + C 2 ) 


Find the minimum and maximum distance from the origin to the ellipse. 


y) — 5z 2 + 6xy + hy 2 — 8 = 0 (A. 25) 

f(x,y) = x 2 + y 2 (A. 26) 

Lagrange Multiplier Rule: 

2x + 2A(5z + 3 y) = 0 (A. 27) 

2y + 2A(3z + 5y) = 0 (A. 28) 

For convenience replace A by — j and divide by 2: 

(5 - A)z + 3y = 0 (A. 29) 

3z + (5 — A )y = 0 Eigenvalue problem (A. 30) 

Solve for A: A! = 2 A 2 = 8 

Inserting in Equation A. 29 or A. 30 yields: 


id = Xi and y 7 = —x 7 <— Insert in Equation A. 25 

From Equation A. 26 one obtains then the minimum and maximum distances: 

d\ = \fx J + yf = \J\J2 + 1/2 = 1 
d% = \fx I + = \/2 + 2 = 2 


211 




References 


Craig, R. R. Jr, “ Structural Dynamics ”, John Wiley and Sons, New York (1981). 

Greenwood, D. T., “ Principles of Dynamics " , 2nd ed., Prentice Hall, Inc., Engle- 
wood Cliffs, N. J. (1988). 

Lanczos, C., “ The Variational Principles of Mechanics” , The University of Toronto 
Press, Toronto (1949). Reprinted by Dover Pulications. 

Meirovitch, L .,“ Methods of Analytical Dynamics ”, McGraw Hill CO., New York 
(1980). 

Sommerfeld, A., “Mechanics" , Academic Press, New York (1952). 

Thomson, W. T., “Introduction to Space Dynamics " , John Wiley and Sons, New 
York (1961). Reprinted by Dover Publications. 


213 


WlgHTHMAICT KAKff 


PRECEDING PAGE BLANK NOT FILMED 



NASA 

National Aeronautics and 
Space Administration 

Report Documentation Page 


1. Report No. 

NASA RP-1262 

2. Government Accession No. 

3. Recipient's Catalog No. 

4. Title and Subtitle 

Methods of Applied Dynamics 

5. Report Date 

May 1991 

6. Performing Organization Code 

7. Author(s) 


8. Performing Organization Report No. 


M. H. Rheinfurth and H. B. Wilson* 

10. Work Unit No 


9. Performing Organization Name and Address 

M-659 

George C. Marshall Space Flight Center 
Marshall Space Flight Center, Alabama 35812 

1 1 . Contract or Grant No. 


13. Type of Report and Period Covered 

12. Sponsoring Agency Name and Address 

Reference Publication 

National Aeronautics and Space Administration 

14. Sponsoring Agency Code 

Washington, DC 20546 



15. Supplementary Notes 


Prepared by Structures and Dynamics Laboratory, Science and Engineering Directorate 
♦Professor of Engineering Mechanics, University of Alabama, Tuscaloosa, Alabama 35847. 

1 6. Abstract ™ ‘ 

The monograph was prepared to give the practicing engineer a clear understanding of dy- 
namics with special consideration given to the dynamic analysis of aerospace systems. It is 
conceived to be both a desk-top reference and a refresher for aerospace engineers in government 
and industry. It could also be used as a supplement to standard texts for in-house training courses 
on the subject. 


17. Key Words (Suggested by Authors)) " 1 18 Distribution Statement 

Kinematics, Euler Angles, Quaternions, Particle 

Dynamics, Rigid Body Dynamics, Gyroscopic Unclassified Unlimited 

Systems, Langrangian Methods, Quasi- 

Coordinates, Modal Synthesis Subject Category 31 



For sale by the National Technical Information Service, Springfield, VA 22161-2171 


NASA-Langley, 1991 











